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1 .  INTRODUCTION 


RPSL1D  is  a  spin  off  of  the  FORTRAN  IV  computer  program  REPSIL1’2 
developed  at  the  Ballistic  Research  Laboratory  to  treat  large,  transient, 
elastoplastic  deformations  in  thin  shells.  Thin  shell  programs,  such  as 
REPSIL,  replace  the  essentially  three-dimensional  geometry  of  the  deform¬ 
ing  shell  with  a  two-dimensional  reference  surface  and  some  assumptions 
about  deformation  through  the  thickness.  The  REPSIL  family  of  programs 
is  based  on  theory  that  restricts  them  to  thin  Kirchoff  shells  with 
negligible  rotary  inertia.  The  coding  restricts  the  programs  to  shells 
of  uniform  thickness  with  clamped  or  hinged  edges  and  a  plane  of  sym¬ 
metry.  The  RPSL1D  version  further  restricts  the  program  to  geometries 
where  the  reference  surface  of  the  shell  is  a  function  of  one  variable 
and  a  symmetry  assumption. 

Three  geometries  have  been  incorporated:  axially  symmetric  shells, 
slab  symmetric  shells,  and  laterally  symmetric  beams.  The  standard 
REPSIL  (with  some  modifications)  could  be  used  for  these  problems,  but 
RPSL1D  is  much  faster  and  uses  about  one-fifth  as  much  memory.  Further, 
RPSL1D  results  should  be  more  accurate;  the  numerical  derivatives  with 
respect  to  the  symmetric  variable  are  replaced  by  exact  derivatives, 
and  spacing  for  the  other  variable  can  be  made  smaller  without  unreason¬ 
able  demands  on  machine  storage  or  time. 

The  simpler  one-dimensional  model  makes  it  easier  to  experiment 
with  alternate  methods  of  computing.  Many  changes  and  additions  have 
been  tried.  Some  of  these  have  been  incorporated  into  RPSL1D,  as  tabu¬ 
lated  in  Appendix  C,  and  others  discussed  in  Section  5  and  Appendices 
E  and  F  may  be  useful  in  the  future. 

The  purpose  of  this  report  is  to  record  RPSL1D,  as  cataloged  in  the 
CDC  CYBER  system  at  BRL,  to  record  some  user  subroutines  and  useful 
options,  to  serve  as  a  user's  guide,  and  to  give  the  basic  information 
for  further  modifications.  This  is  not  a  complete  guide;  the  user 
should  also  obtain  a  copy  of  Reference  2,  the  user's  manual  for  REPSIL. 
An  understanding  of  the  theory  given  in  Reference  1  is  also  very  help¬ 
ful.  The  symbols  in  the  present  report,  which  are  sometimes  used  with¬ 
out  sufficient  explanation,  are  those  of  Reference  1. 

Sections  2  and  3  and  the  list  of  equations  in  Appendix  A  should 
bridge  the  gap  between  REPSIL  and  RPSL1D.  In  an  effort  to  keep  this 
report  within  reasonable  bounds,  we  will  not  repeat  theory  given  in 
Reference  1  nor  will  we  repeat  much  of  the  information  given  in  the 
REPSIL  user's  manual.  Reference  2. 


M.  Santiago ,  "Formulation  of  the  Large  Deflection  Shell  Equations 
for  Use  in  Finite  Difference  Structural  Response  Computer  Codes, " 

U.S.  Army  Ballistic  Research  Laboratories,  Report  No.  1571,  Feb.  72. 

2j.  M.  Santiago,  H.  L.  Wisniewski,  N.  J.  Buffington,  Jr.,  "A  User's 
Manual  for  the  REPSIL  Code,"  U.S.  Army  Ballistic  Research  Laboratories, 
Report  No.  1744,  Oct.  74. 
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The  RPSL1D  user  must  prepare  input  as  described  in  Section  6  and 
must  supply  proper  "user"  subroutines  as  described  in  Section  4.  He 
may  retain  the  user  subroutines  with  the  cataloged  program  listed  in 
Appendix  C,  replace  them  with  copies  of  user  subroutines  listed  in 
Appendices  E  and  F,  or  replace  them  with  his  own  subroutines. 

For  some  particular  applications  it  seems  inevitable  that  some 
modifications  are  needed  or  desirable.  A  number  of  possibly  useful 
options  are  described  in  Section  5  and  listed  in  Appendices  E  and  F. 

The  program,  as  stored  in  UPDATE  form  in  the  CDC  7600,  is  listed 
in  Appendix  C.  A  complete  list  of  FORTRAN  names  used  in  RPSL1D  is  given 
in  Appendix  B  along  with  their  definition  and  some  explanation.  Two 
example  problems  are  discussed  in  Appendix  E. 

The  accessory  program  for  plotting  output  from  RPSL1D  is  discussed 
and  listed  in  Appendix  D. 

The  main  purpose  of  this  report  is  to  serve  as  a  guide  to  the  use 
and  possible  modification  of  the  RPSL1D  program.  The  author's  interest 
in  this  program  is  as  a  programmer.  Hence,  there  is  a  bias  in  the  pre¬ 
sentation  toward  the  frequent  use  of  FORTRAN  names  and  formulation,  and 
the  trivial  details  of  programming,  at  the  expense  of  physical  meaning. 


2 .  MODELS 

The  first  one-dimensional  version  of  REPSIL  treated  axially  sym¬ 
metric  shells.  This  was  modified  to  handle  slab  symmetric  shells,  and 
this  in  turn  changed  to  treat  symmetric  beams. 

The  position,  x*,  of  any  point  in  the  shell  may  be  expressed  in 
terms  of  a  point,  £,  on  the  reference  surface,  the  normal,  n,  to  the 
reference  surface  at  £,  and  the  distance  from  the  reference  surface, 

x  =  r  +  cn  .  (1) 

Points  on  the  reference  surface  are  expressed  in  Cartesian  coordinates 
as 

I  =  yJ  ij  5  1,1  ii  ♦  y2  *2  *  y5  ij  • 

where  the  are  the  orthonormal  basis  for  vectors  in  3-space.  The  , 

and  hence  other  dependent  parameters,  are  functions  of  the  material 
(Lagrangian)  parameters,  and  F,2 ,  and  of  time,  t. 

r  =  r(^}  ,£2,t)  =  (51  ,£;2,t)  i^  .  (3) 


The  symbols  of  Reference  1  will  usually  be  followed  in  this  report. 
The  tilde  helots  a  letter  denotes  a  vector  (e.g.  gc). 


Our  characterization  of  axial  and  slab  symmetry  is  through  speciali¬ 
zation  of  the  vector  x.  and  the  material  coordinates  C1  and  C2.  The  basic 

specializations  are  given  in  this  section.  A  summary  of  equations  is 
given  in  Appendix  A. 

The  type  of  model  is  prescribed  for  RPSL1D  by  control  parameters 
assigned  in  the  user's  subroutine  INGEOM. 

2.1  Axially  Symmetric  Shells 

An  axially  symmetric  shell  is  a  shell  of  revolution  (see  Figure  2.1) 
such  as  a  circular  cylinder  or  the  surface  of  a  truncated  right  circular 
cone.  A  shell  of  revolution  is  the  surface  generated  by  revolving  a 
plane  curve  (actually  two  parallel  curves  so  the  shell  has  thickness) 
about  a  line,  called  the  axis  of  revolution,  in  its  plane.  Every  point 
of  the  revolving  curve  describes  a  circle  whose  center  is  on  the  axis. 

For  our  shell,  the  curve  must  be  smooth  and  must  not  intersect  the  axis 
in  the  region  of  interest. 

We  study  the  motion  of  the  shell  by  following  the  motion  of  the 
curve  formed  by  the  intersection  of  the  shell's  reference  surface  and 
a  plane  through  the  axis.  We  do  not  allow  twisting  in  our  model,  so 
the  initial  motion,  and  any  loading,  must  be  in  this  plane.  Because  of 
the  axial  symmetry,  the  internal  forces  are  all  in  this  plane. 

Let  51  =  0,  C2  =  C,  and  introduce  the  orthonormal  base  vectors 
(£>.£»£)  which  vary  with  angle  0.  These  are  related  to  (ij,^,^)  bY 
the  relations 


r  =  r(0)  =  cos0  +  sin0  , 

(4) 

£  =  0,(0)  =  -sin0  i3  +  cos©  ij  , 

(5) 

i=  k2  ■ 

(6) 

Then, 

X  *  £(e,5,t)  =  R  i  +  z  i  (7) 

where 

R  =  RU,t)  =  {(Y3)2  +  (Y1)2}*5  ,  (8) 

Z  s  ZU,t)  =  Y2  .  (9) 
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Figure  2.1  -  Model  for  axial  symmetry  showing  the  two 
orthonormal  bases  (ij.i^.ij)  and  (£,&,]<)  with  origin  at 
the  center  of  end  1  of  the  shell,  and  a  vector  £  to  the 
surface  of  the  axisymmetric  shell. 

The  same  formulation  holds  for  every  angle  6.  Since  evaluation 
is  needed  at  only  one  angle,  we  naturally  choose  0=0.  Then, 


To  test  this  new  program  for  axial  symmetry,  the  same  problem  (an 
initially  cylindrical  shell  subjected  to  a  uniform  impulsive  load)  was 
simulated  with  both  programs.  Both  programs  were  run  undamped  for  200 
time  cycles  and  then  damped.  Both  self  terminated  on  time  cycle  262. 
The  results  were  very  close.  The  maximum  difference  noted  in  displace¬ 
ment  was  less  than  2*10~8.  Since  the  displacement  was  about  unity, 

2*10  8  was  also  the  relative  difference  in  displacement.  Strains  and 
displacement  increments  are  ,ery  sensitive  to  program  differences.  The 
largest  relative  difference  noticed  in  either  was  10"5.  (There  were 
relative  differences  of  about  0.003  in  all  the  printed  energy  values. 
This  was  because  the  circumference  of  the  right  circular  cylinder  was 
a  circle  in  the  axial  symmetric  program  but  was  a  48-sided  regular 
polygon  in  the  standard  REPSIL.) 

2.2  Slab  Symmetric  Shells 

Slab  symmetric  shells  may  be  characterized  by  the  displacement  of 
any  point  being  a  function  of  its  initial  distance,  52,  from  end  1  (see 
Figure  2.2). 

Let  Y1  =  5l,  Y2  =  Y2U2,t),  and  Y3  =  Y3U2,t).  Then, 

r(Cl,S2,t)  =  +  Y2(C2,t)  i2  +  Y3(52,t)  i3  .  (11) 

To  utilize  the  formulation  and  programming  for  axial  symmetry,  we 
renamed  variables 

S2  =  y3  =  R,  and  Y2  -  Z  ,  (12) 

and  confined  evaluation  to  the  plane  Y1  =  0.  Then, 

X  =  X(0,5,t)  =  R(€,t)  i3  +  Z(5,t)  i2 

=  Y3(?,t)  i3  +  Y2(C,t)  i2  .  (13) 


Figure  2.2  -  Model  for  Slab  Symmetry. 


To  maintain  slab  symmetry,  the  loading  must  not  contain  any  compo¬ 
nent  in  the  ^  direction. 

For  both  axial  and  slab  symmetry,  2  is  a  distance  along  the  ^ 

axis.  For  axial  symmetry,  R  is  the  radial  distance  from  the  axis  of 
rotation.  For  slab  symmetry,  R  is  the  perpendicular  distance  from  the 

=  0  plane. 

The  initial  slab  symmetric  coding  was  also  tested  by  comparing  with 
results  from  a  standard  REPSIL  run.  (A  minor  coding  change  was  needed 
to  insert  cyclic  boundary  conditions  in  REPSIL.)  The  results  were 
nearly  identical.  The  maximum  relative  errors  in  any  variables  were  of 
the  order  of  10  13 . 


2.3  Laterally  Symmetric  Beams 

One  principal  difference  between  slab  symmetric  shells  and  beams 
is  in  the  change  of  concept  from  area  to  unit  length  (e.g.  pressure, 
i.e.  force/unit  area,  for  shells  and  force/unit  length  for  beams.)  The 
original  beam  considered  was  a  narrow  rectangular  beam.  It  was  treated 
as  a  special  case  of  slab  symmetry  with  a  uniaxial  stress-strain  relation 
replacing  the  biaxial  one.  The  uniaxial  stress-strain  was  imposed  by 
using  a  new  subroutine  BMSTRS,  instead  of  subroutine  STRESS,  for  beams. 

By  setting  the  width,  DETA1  (A^1),  to  unity,  the  concept  of  surface  area 
became  length  without  any  further  changes. 

Treating  beams  that  were  not  rectangular  was  more  involved.  It 
required  changes  in  fact  as  well  as  in  concept.  The  main  change  was 
for  integrals  through  the  thickness.  The  simplest  form  of  numerical 
integration  through  the  thickness  is  used  in  the  standard  REPSIL: 


f(0  dc 


The  thickness,  h,  is  divided  into  L  layers  of  equal  thickness,  A£  =  h/L, 
and  the  function  to  be  integrated  is  evaluated  at  each  where 

is  in  the  center  of  the  k'th  layer.  The  parameter  x,  is  measured  from 
the  reference  surface  in  the  center  of  the  shell. 

For  a  beam,  an  integral  through  the  thickness  is  replaced  by  an 
integral  over  the  cross-sectional  area.  When  multiplied  by  the  proper 
metric,  this  becomes  the  quantity  per  unit  surface  area.  For  a  beam, 
the  desired  result  is  the  quantity  per  unit  length.  For  a  rectangular 
beam  we  could  simply  multiply  the  previous  result  by  the  width  of  the 
beam,  A?1: 

//•+h/2  L 

f(c)  dA  =  A?1  /  f(0  d?  ~  A?1  A?  £  f(Ck)  •  (15) 

A  J -h/2  k=l 
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This  has  been  replaced  by  a  more  general  form  for  laterally  sym¬ 
metric  but  not  necessarily  rectangular  beams: 

/  f(0  “A  '  g,  \  f‘V  (16) 

where  the  ^  are  not  necessarily  evenly  spaced  and  is  the  area  assoc¬ 
iated  with  The  ?k  and  should  be  assigned  so  that  is  the 

cross-sectional  area  of  the  beam,  Wk  =  0  (so  the  reference  surface 

is  the  neutral  surface  for  bending),  and  £(Ck)2  Wk  is  the  area  moment 

of  inertia  of  the  cross  section.  (This  form  of  integral  approximation 
is  also  suitable  for  Gaussian  integration,  which  is  preferred  for  shells; 
see  Section  3.4.) 

See  the  discussion  of  INGEOM  for  a  beam  (Section  4.1.1)  for  a 
particular  example.  The  user  subroutine  INGEOM  for  beams  must  also 
supply  the  upper  limit  of  the  beam,  ZU  =  cu.  Subroutine  START  assigns 
=  ZL  =  ZU  -  h.  For  beams,  ZU  and  ZL  merely  locate  positions,  upper 

and  lower  respectively,  at  which  surface  strains  are  computed  for  out¬ 
put.  The  shell  versions  set  ZU  =  h/2,  ZL  =  -h/2,  and  ZU  is  used  in 
START  to  compute  a  number  of  program  constants. 

Results  from  RPSL1D  for  small  elastic  displacements  for  vibrations 
initiated  in  the  fundamental  transverse  modes  for  clamped,  hinged,  and 
cantilever  beams  are  in  good  agreement  with  simple  beam  theory.  An 
example  of  buckling  from  lateral  forces  on  the  ends  also  agreed  with 
theory . 


3.  CHANGES  FROM  REPSIL 

The  principal  reason  for  programming  a  one-dimensional  version  of 
REPSIL  was  to  have  a  simple  vehicle  for  testing  possible  changes  before 
implementing  them  in  REPSIL.  However,  the  program  proved  very  useful 
for  a  variety  of  problems.  Many  changes  have  been  coded.  Some  of 
these  were  purely  experimental,  but  most  of  them  were  for  particular 
problems.  Most  of  the  changes  have  been  superseded  or  were  not  con¬ 
sidered  general  enough  to  transfer  to  the  CDC  computer.  Some  of  the 
changes  have  been  incorporated  as  part  of  RPSL1D,  as  listed  in  Appendix  C, 
and  some  are  retained  as  options. 

In  this  section  we  will  outline  the  changes  made  to  REPSIL  to 
create  the  RPSL1D  program  stored  in  the  CDC  7600  as  of  September  1978, 
and  describe  the  new  boundary  condition  implementation,  the  new  dif¬ 
ferencing  scheme,  Gaussian  integration,  and  the  new  stability  criteria. 


3.1  Outline  of  the  Development  of  RPSL1D 

Initially,  the  REPSIL  program,  listed  in  Reference  2,  was  changed  to  a 
one-dimensional  code  for  axially  symmetric  shells  (see  Section  2.1). 

A  change  of  notation  was  introduced  to  eliminate  errors. 

This  program  was  modified  to  include  slab  symmetry  (Section  2.2). 

This  was  then  modified  to  include  beams  with  a  vertical  axis  of 
symmetry  (Section  2.3). 

The  method  of  imposing  boundary  conditions  was  changed  by  inserting 
artificial,  external  points  as  described  in  Section  3.3,  and  the 
divided  difference  approximations  for  partial  derivatives  were  changed 
(see  Section  3.2)  to  a  form  using  only  central  differences.  The 
system  was  made  more  compact  by  evaluating  at  mid-mesh  points. 

The  code  was  changed  to  permit  symmetry  at  either  end  and  permit  a  free 
end  at  end  2.  This  change  was  intended  for  beams  only.  (An  option 
called  APLFRC,  see  Section  5.4,  which  permits  forces  to  be  applied 
at  either  or  both  ends  of  beams  is  available.) 

Coding  was  added  to  monitor  maximum  deflection  and  extreme  surface 
strains,  and  to  print  these  extremes  along  with  the  surface  strain 
prints. 

The  computation  of  surface  strains  was  revised.  The  original  surface 
strain  computation  was  a  copy  of  that  described  on  pages  50  through  53 
of  Reference  2  except  that  a12  and  b12  were  zero  and  ^  and  ^  replace 

h/2  and  -h/2.  This  was  changed  so  that  the  interpolation  in  the  non¬ 
zero  covariant  components  of  strain,  and  e22  (e^  =  ^(EAa^  ±  h  EAb^)) 

was  replaced  by  interpolation  in  the  changes  in  the  covariant  components 
of  the  middle  surface  metric,  EAa^  and  EAa22,  and  of  the  second  funda¬ 
mental  tensor,  EAb^  and  IAb22.  Component  EAa22  is  computed  between 
mesh  points . 

Optional  Gaussian  integration  through  the  thickness  of  shells  was 
introduced . 

Several  minor,  mostly  cosmetic,  changes  were  made. 

The  coding  was  transferred  from  the  BRLESC  computer  to  the  CDC  7600 
(with  more  minor  changes)  and  stored  in  the  UPDATE  form. 

3.2  New  Difference  Operators 

The  primary  purpose  of  the  REPSIL  family  of  programs  is  to  solve  the 
equations  of  motion.  This  is  done  by  a  sequence  of  equations  (see  Reference  1, 
pp.  99-106)  that  carry  the  solution  forward  one  time  step. 
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Difference  approximations  of  derivatives  with  respect  to  the  material  coor¬ 
dinates,  C1  and  C2,  are  required  at  two  points  in  this  sequence  of  equations. 
First,  we  need  approximations  to  the  first  and  second  derivatives  of  position, 

£,  and  the  incremental  displacement,  A£,  to  define  the  local  geometry.  Finally, 

we  need  approximations  for  first  partial  derivatives  of  the  components  of  the 
stress  resultant  tensor,  U*a,  and  second  partial  derivatives  of  the  normal 

~  ot  ft 

vector  components  of  the  bending  resultant  tensor,  J4*  ,  for  the  equations  of 

motion. 

The  REPSIL  equations  of  motion  for  a  mesh  point,  as  given  in  Reference  1, 


.  +  .  (At)2  r  32bi  8  3N*  . 

Aji  =  Aji  +  +  £  .  (17) 

Y*o  aeV8  3Ca 

The  change  in  the  displacement  of  the  point  from  time  t  to  t  +  At  is  Ajj+. 

The  position  at  time  t  +  At  is  £+  =  £  +  Aj^+,  where  £  is  the  position  at  time  t. 

The  mass  per  unit  initial  middle  surface  area,  y*^,  is  constant  for  a  given 

point.  (The  term  in  brackets  has  units  of  force  per  initial  unit  area.) 
v 

P*  *  -P  a  n  is  the  force  per  initial  unit  area  along  the  normal,  n,  due  to  a 

pressure  P.  (a  is  the  determinant  of  the  covariant  metric  tensor.)  The  one¬ 
dimensional  assumptions  remove  a  number  of  the  terms  in  the  equations  of 
motion.  There  are  no  terms  in  the  i^  direction,  M*12  =  M*21  =  0,  and  the 

derivatives  with  respect  to  are  simple  functions  which  are  computed 
explicitly.  The  equations  of  motion  become 


A  +  .  (At)2  r  32M*22  3N*2  . 

ax  =  a£  +  2  +  £ 

Y*o  I  oc2)  n2 

-CM*11  *  i2  +  N*1  .  ij)  i2l 
+  .  +  J 


r  *  t  ♦  •  (19) 

(The  transformation  to  this  equation  is  outlined  in  Appendix  A.)  The  vectors 
JSj*22,  ij*2,  and  their  derivatives  and  £*  and  may  all  have  components  in 

both  the  i3  and  £2  directions.  JJ*1  has  a  component  in  the  direction  only. 

For  slab  symmetry  and  beams,  the  term  (Hi**1  •  i2  +  K**  *  ij)  *s  deleted. 

For  beams,  y*Q  is  mass/unit  length  and  the  term  in  brackets  is  force/unit  length. 


The  difference  scheme  used  in  the  standard  REPSIL  has  two  defects  which 
have  occasionally  been  troublesome:  the  use  of  forward,  or  backward,  dif¬ 
ferences  at  fixed  edges  may  introduce  a  slowly  developing  instability  which 
cannot  be  eliminated  by  a  reduction  in  the  time  step.  (The  presence  of  this 
instability  could  be  detected  by  an  eigenvalue  analysis  of  the  linearized 
equations  of  motion.  We  have  carried  out  this  analysis  for  a  few  cases  with 
very  small  grids.  Unfortunately,  for  most  useful  meshes,  the  eigenvalue 
analysis  would  take  longer  than  the  REPSIL  program.)  The  other  defect  is 
that  the  differences  for  first  derivatives  are  not  compact  enough.  (As  an 
experiment  with  the  initial  RPSL1D  code,  alternate  points  in  a  slab  symmetric 
run  were  at  rest  or  given  an  initial  axial  velocity.  The  points  initially 
at  rest  remained  at  rest.  The  motion  of  the  moving  points  produced  a  reason¬ 
able  solution.)  These  two  defects  have  been  removed  from  RPSL1D  (and  from 
a  modified  version  of  REPSIL)  by  modeling  fixed  boundaries  with  artificial 
external  points  and  by  introducing  new  difference  operators  which  compute 
and  use  forces  between  mesh  points  as  well  as  at  mesh  points. 

Three  sets  of  central  difference  operators  are  used  in  RPSL1D  to 
compute  first  and  second  differences  with  respect  to  the  Lagrangian  coor¬ 
dinate  £  =C2.  To  simplify  the  notation,  replace  £2  by  y  and  use  the  notation 
fy(n)  and  fyy(n)  to  represent  the  approximations  to  3f/3£2  and  32f/(9£;2)2. 

respectively,  evaluated  at  the  n'th  mesh  point. 

The  first  set  of  difference  operators  is  used  in  subroutine  GRAD  to 
approximate  derivatives  of  position  components,  and  derivatives  of  compo¬ 
nents  of  displacement  increments,  at  mesh  points: 

fy(n)  =  (f(n+l)  -  f(n-l)}/(2  Ay)  ,  (2°) 

f  (n)  =  (f (n+1)  -  2  f (n)  +  f (n-1) }/(Ay) 2  .  (21) 

The  second  set  of  difference  operators  is  used  in  subroutine  GRAD  to 


approximate  the  same  derivatives  at  mesh  midpoints: 

fy(n+%)  =  (f (n+1)  -  f (n) }/ Ay  ,  (22) 

fyyCn+%)  =  %{f (n+2)  -  f (n+1)  -  f(n)  +  f (n-1) }/(Ay) 2  .  (23) 

The  third  set  of  central  difference  operators  is  used  in  subroutine  MOTION 
to  compute  terms  for  the  equations  of  motion: 

fy(n)  =  {f(n+y  -  f  (n-4) )/ Ay  ,  (24) 

f  (n)  =  {f (n+1)  =  2  f (n)  +  f (n-1) }/(Ay) 2  .  (25) 


These  difference  operators,  which  are  used  at  fixed  edges  as  well  as  at 
internal  points,  are  all  central  difference  operators.  The  resulting  system 
of  difference  equations  is  stable  with  the  proper  choice  of  time  step.  At. 
Unfortunately,  when  the  program  is  run  for  a  beam  with  a  free  end  at  end  2, 
say  n  =  N,  backward  difference  operators  have  to  be  used  at  the  free  end. 

In  GRAD,  for  a  free  end,  we  use: 
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f  (N)  =  {3  f (N)  -  4  f(N-l)  +  f (N-2) }/ (2  Ay)  ,  (26) 

£  (N)  =  (2  f (N)  -  5  f(N-l)  +  4  f (N-2)  -  f (N-3) }/(Ay) 2  ,  (27) 

fy(N-y  =  (f (N)  -  f (N-l) }/Ay  ,  (as  before)  (28) 

fyyfN-Js)  =  55(3  f  (N)  -  7  f  (N-l)  +  5  f(N-2)  -  f  (N-3)  }/ (Ay)  2  .  (29) 

In  MOTION,  for  a  free  end,  we  have: 

f  (N)  =  {f (N)  -  f(N-3s)}/(Ay/2)  ,  (30) 

f  (N)  =  (2  f (N)  -  5  f (N-l)  +  4  f (N-2)  -  f (N-3) }/ (Ay) 2  .  (31) 

It  is  also  necessary  to  approximate  N*2  =  -SM*22/3£  at  the  free  end.  This 
is  done  in  subroutine  RESULT  with 

N*2(N)  =  Jj{4  M*22(N-1)  -  M*22 (N-2)  }/ A?  ,  (32) 

because 

M*22(N)  =  0.0  .  (33) 

3.3  External  Points  and  End  Conditions 


The  end  conditions  for  a  fixed  or  symmetric  end  is  imposed  in  RFSL1D 
by  the  initial  positioning  and  the  incrementing  of  a  fictional  external 
point.  The  following  notation  is  used  to  explain  this: 

£  A  vector  to  a  point  on  the  reference  surface  at  time  t. 

(£  -  R  i3  +  Z  i2) 

r+  A  vector  to  a  point  on  the  reference  surface  at  time  t  +  At. 

*X  X+  -  X 

Tg  A  vector  to  the  end  point. 

Tj  A  vector  to  the  internal  mesh  point  adjacent  to  r^. 

r£  A  vector  to  the  external  point  adjacent  to  £g. 

ng  The  unit  normal  to  the  reference  surface  at  £g. 

(nD  =  n  i_  +  n,  i_) 

*-B  r  ~3  k  ~2 

i  The  unit  tangent  to  the  reference  surface  at  rfi. 

<i  -  -nk  A3  *  nr  i2) 

The  initial  position  of  external  points  is  established  in  subroutine 
BOUNDR  and  the  corresponding  increments  are  inserted  in  BOUNDU.  The 
initial  normal  at  clamped  end  points  is  supplied  in  a  new  user  subroutine 
INNORM. 


3.3.1  Hinged  Ends, 
point  at  a  hinged  end  is 


The  increment  is 

ArE  -  -toj  .  (35) 

We  also  set 

Arg  =  0  .  (36) 

These  inserts  impose  r^  is  constant  and  M*22  =  0  at  the  boundary.  These 
are  the  conditions  for  a  hinged  end  as  given  in  Reference  1. 

3.3.2  Clamped  Ends.  The  initial  position  of  an  artificial  external 
point  adjacent  to  a  clamped  end  is 

r_  =  rT  -  2{(r.  -  rD)  •  t)  t  .  (37) 

The  corresponding  increment  is 

aIe  “  ATj  -  2{Arj  •  T)  t  .  (38) 


The  initial  position  of  an  artificial  external 


~E  =  2  ~B  "  ~I  ' 


(34) 


We  also  impose 


A£B  =  0  . 


(39) 


These  inserts  force  rg  and  ng  to  remain  constant.  These  are  the  condi¬ 
tions  for  a  clamped  end  given  in  Reference  1.  (The  components  of  ng  are  set 
in  a  user  subroutine  INNORM  and  never  changed  for  a  clamped  end.) 


3.3.3  Symmetric  Ends.  The  treatment  of  symmetric  ends  is  simply  a 
reduction  of  that  for  symmetric  edges  in  REPSIL.  A  symmetry  plane  is 
assumed  perpendicular  to  the  basis  vector  .i^  The  initial  position  of  an 

external  point  adjacent  to  a  symmetry  end  is 

%  =  ~I  2^£i  ■  *  ^2*  ~2  *  (40) 

The  change  in  r^  at  each  time  step  is 

A£e  «  A£j  -  2{ArI  •  i2)  i2  .  (41) 


-  =  T. 
2  ~ 


These  are  the  same  equations  used  for  clamped  ends  with  i 
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The  external  points  at  symmetry  ends  represent  real  points.  They  are 
not  artificial  like  those  outside  the  fixed  ends.  Just  as  in  the  standard 
REPSIL,  we  need  to  know  membrane  and  bending  resultant  components  (N*2  and 
M*22,  respectively)  outside  the  symmetry  end.  We  have 


(M*22)j  =  M*22  nj  =  (FM22) j  i3  +  (FM22)j  i2  ,  (42) 

defined  at  rT,  and 

=  (fnJ)j  ix ,  (43) 

(N*2)j  =  (FN2)j  i3  +  (FN2)j  i2  ,  (44) 

defined  between  rT  and  rD.  The  external  values  needed  are 
~I 

(FMj2)E  =  (FM22) J  ,  (45) 

(FM22)e  *  -(FM22)j  ,  (46) 

(FNg)E  =  CFNq)i  .  (47) 

(FNj)E  =  -(FN2)j  ,  (48) 

(FN2)e  =  (FN2)  j  .  (49) 


3.3.4  Free  Ends.  One  of  the  first  problems  posed  for  the  RPSL1D  code 
was  for  a  cantilever  I-beam.  To  model  this  beam  a  free  end  was  needed.  The 
conditions  for  a  free  end  (see  Ref.  2,  Eq.  6.52)  are  imposed  by 


M*22  =  0.0  ,  N*2  =  - 9M*22/A£  . 


(50) 


Regrettably,  it  is  necessary  to  use  non-central  differences  to  approximate 
this  derivative,  and  other  derivatives,  at  free  ends  (see  Section  3.2) . 

The  cataloged  program  permits  end  2  to  be  free  but  not  end  1.  (The 
optional  coding  called  APLFRC,  see  Section  5.4,  allows  either  end  to  be 
free  with,  or  without,  applied  forces.  The  applied  forces,  or  absence 
of  forces,  are  supplied  by  a  user’s  subroutine  called  ENDFRC.) 

3.4  Optional  Gaussian  Integration 

The  original  RPSL1D  programming  approximated  integration  through  the 
thickness  of  shells  by  Riemann  sums  assuming  equal  layers  about  evenly 
spaced  points.  This  is  replaced  with  Gaussian  integration  only  if  IGAUSS 
is  1  and  IB  is  not  1.  (These  control  parameters  are  set  in  the  user  sub¬ 
routine  INGEOM,  see  Section  4.1.)  IB  =  1  specifies  a  beam.  We  do  not  use 
Gaussian  integration  with  beams.  We  use  Riemann  sums  with  unevenly  spaced 
points  and  "weights"  equal  to  their  associated  cross-sectional  area  (as 
explained  in  Section  2.3).  For  the  Riemann  sums  with  L  equal  layers  we 
choose 
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ck  =  h  h  (1  -  (2  k  -  1)  /  L)  ,  (k  =  1,2, ---.L) 

AC  =  h/L  , 

Si  =  h/2  * 

=  -  h/2  , 


and  approximate  integration  through  the  thickness  by 

f^U  f (?)  dc  ~  £  f(cJ  AC., 

•'c.  k=l 


where  h  is  the  shell  thickness,  ?  and  ?  are  the  distances  from  the  reference 

surface  (the  middle  surface)  to  the  upper  and  lower  surfaces,  respectively. 

L  is  the  number  of  equally  spaced  layers,  and  ?k  is  the  center  of  the  k'th 
layer. 

For  Gaussian  integration,  we  replace  the  ?k  by 

?k  “  ^  h  ,  (56) 


Wk  "  WkL  > 


and  use  the  approximation 


Jr?  K 

[  u  f(c)  dc  ~  (h/2)  £  f(ck)  Wk  . 

k=l 


where  the  x^^  and  are  from  tables  for  Gaussian  integration  (e.g.  Table 

25.4,  Handbook  of  Mathematical  Functions,  National  Bureau  of  Standards).  The 
Xj^  are  zeros  of  the  Legendre  Polynomial 

p  (x,  .  1  ^  Cx2  -  1J.1-  (=9) 

L  2L  L!  (dx)L 

Two  sets  of  integrals  are  affected:  the  integral  for  strain  energy  in 
subroutine  STRESS,  and  the  integrals  for  force  and  moment  resultants  in 
subroutine  RESULT. 

The  Riemann  sums  for  the  force  and  moment  resultants  at  a  point  N  are 
of  the  form 

IN  =  TA  £  fN(ck)  ,  (60) 

1c 
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TA  =  A?  =  h/L  . 


(61) 


The  signifies  the  value  of  the  appropriate  function  at  mesh  point  N 

and  integration  station  where  the  ^  are  in  the  center  of  evenly  spaced 
layers. 

All  that  is  needed  to  change  these  sums  to  Gaussian  integration  is  to 
select  the  proper  and  W^,  replace  fN(?k)  by  Wk  £N^k^’  and  set  TA  to  h^2. 

For  example. 


h  - (h/2)  i  \  w  • 

k=l 


(62) 


(The  form  is  the  same  for  approximations  of  integrals  over  cross-sectional 
area  for  beams.  This  is  done  with  Riemann  sums  but  with  unevenly  spaced 

The  and  equal  to  the  area  associated  with  ,  are  selected  in  the 

user  subroutine  INGEOM,  and  TA  is  set  to  1.0  in  START.) 

3.5  Convergence  and  Stability 

Some  examination  of  convergence  and  stability  is  essential  for  any 
numerical  solution  of  differential  equations.  The  usual  method  of  examining 
convergence  is  to  find  solutions  for  a  particular  problem  using  an  increasing 
number  of  meshes  until  these  solutions  converge.  Such  tests  with  the  REPSIL 
family  of  programs  have  shown  remarkable  insensitivity  to  the  number  of 
meshes  for  problems  with  large  plastic  deformation.  In  those  cases  tested 
it  seems  that  the  maximum  deflection,  usually  at  the  center,  decreases 
slightly  with  additional  meshes,  and  the  strain  and  deflection  increase 
slightly  near  the  edges.  For  elastic  vibration  problems,  more  meshes  may 
be  needed  to  permit  high  vibration  modes.  The  choice  of  mesh,  and  hence 
the  accuracy,  must  ultimately  be  decided  by  the  user. 

Numerical  solutions  of  differential  equations  by  explicit  methods  such 
as  that  used  with  the  REPSIL  programs  will  be  unstable,  because  any  intro¬ 
duced  error  grows,  unless  the  increment  in  the  independent  variable.  At  in 
our  case,  is  small  enough.  Considerable  time  and  effort  have  been  expended 
in  studying  the  stability  for  the  various  REPSIL  codes.  Even  a  cursory 
discussion  of  the  stability  analysis  is  beyond  the  scope  of  this  report. 

The  following  bounds,  based  on  a  linear  analysis  of  the  non-linear  difference 
equations  of  motion  for  RPSL1D  assuming  the  central  difference  operators, 
are  used  in  RPSL1D: 


Atj^  =  AC/C  ,  (63) 

At-  =  ^s(AC)  (CC~)  .  (64) 


Here,  At^  and  Atg  are  the  stable  time  step  for  the  linearized  membrane  and 
bending  equations  of  motion,  respectively.  C  is  the  longitudinal  wave  speed. 

C  =  (E/p)  2  ,  for  beams.  (65) 

C  =  {E/  p(l  -  v2) ,  for  shells.  (66) 


Here,  E  is  Young’s  modulus  of  Elasticity,  p  is  mass  density, 
Poisson's  ratio.  is  the  approximation  of  {  /  £2  dA  /  / 

A  A 


For  beams. 


and  v  is 
dA}*5  . 


c2  -  <  E  \  (sk’2/  £  v’5  ■ 


(67) 


For  shells  h  units  thick,  using  Riemann  sums  with  L  equal  layers, 

C2  =  {(1  -  l/l2)  h2/^}*5  . 


(68) 


For  shells  with  Gaussian  integration  the  approximation  is  exact: 

C2  =  (h2/12)h  . 


(69) 


The  program  chooses  a  two  digit  truncated  value  0.95  times  the  minimum 
of  AtM  and  Atg.  This  is  used  instead  of  the  input  At  on  input  card  3  if  it 

is  smaller  than  the  input  value  given,  or  if  the  input  value  is  negative  or 
zero. 


The  term  AC  used  by  RPSL1D  in  these  equations  is  the  DETA2  supplied  by 
the  user  through  subroutine  INGEOM.  The  equations  produce  critical  values 
for  At^  and  Atg  if  A?  is  the  minimum  distance  between  mesh  points.  Both  At^ 

and  Atg  have  been  shown  to  be  close  to  the  limits  for  stability  for  runs 

with  fixed  ends  and  a  linearly  elastic  stress-strain  relation.  The  user 
should  examine  all  results  for  numerical  instability.  This  shows  up  first 
as  rapid  vibrations  in  strain  plots.  The  inclusion  of  plasticity  tends  to 
damp  these  vibrations,  but  they  can  still  be  detected.  The  user  should  be 
particularly  critical  if  there  is  a  free  or  forced  end  (non-central  differ¬ 
ences  are  used),  or  if  the  structure  is  compressed.  The  stability  formuli 
are  based  on  the  assumption  that  AC  is  the  minimum  distance  between  mesh 
points.  The  factor  0.95  permits  some  compression,  but  only  up  to  5%  com¬ 
pression  if  At^  is  critical,  or  2.5%  if  Atg  is  critical.  We  have  not  been 

able  to  prove  the  stability  of  the  RPSL1D  solutions  when  non-central  differ¬ 
ence  equations  are  used,  but  we  have  not  detected  instability  when  a  time 
step  that  satisfies  the  stability  criteria  is  used. 
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4.  USER  SUBROUTINES 

There  are  four,  or  possibly  five,  subroutines  in  RPSL1D  which  the  user 
must  either  write,  or  select  from  existing  versions.  These  "user"  sub¬ 
routines  describe  the  geometry  and  loading.  Three  of  them,  INGEOM,  PRESS, 
and  INVEL,  are  one-dimensional  versions  of  subroutines  described  in  Refer¬ 
ence  2.  Subroutine  INNORM  was  added  to  supply  the  unit  normal  to  the 
reference  curve  at  clamped  ends.  Subroutine  ENDFRC  must  be  included,  to 
supply  the  forces  and  moments  on  the  ends  of  a  beam,  if  the  optional  coding 
APLFRC,  which  assumes  this  type  of  loading  (see  Section  5.4),  is  added  to 
the  RPSL1D  program.  Some  of  these  subroutines  are  always  called  and  some 
are  called  under  control  of  input  parameters.  Any  of  the  user  supplied 
subprograms  may  require  input  (this  is  under  the  control  of  the  user) .  The 
input  for  them  follows  all  the  other  input  (input  is  discussed  in  Section  6) . 
The  user  subroutines  are  listed  below  in  the  order  they  are  called,  with  the 
conditions  for  calling. 

INGEOM  is  called  from  START  shortly  after  the  start  of  each  run. 

INNORM  is  called  immediately  after  INGEOM  from  START. 

PRESS  is  called  from  the  initiation  part  of  the  main  program, 

RPSL1D,  if  LOAD  #  0.  It  will  be  called  each  time  cycle 
thereafter  as  long  as  LPRESS  >  NCYCLE. 

INVEL  is  called  in  the  initiation  sequence  in  RPSL1D  if  LOAD  <  0. 

ENDFRC  is  called  from  RESULT  each  cycle,  except  cycle  0,  if  the 

APLFRC  optional  coding  is  included  in  the  program. 

Information  is  brought  to  each  of  the  user  subroutines,  and  carried 
back  to  the  rest  of  the  program,  through  COMMON.  The  RPSL1D  program  is 
intended  to  be  stored  in  the  UPDATE  format.  All  the  COMMON  is  contained  in 
a  COMDECK  called  MAIN  which  is  listed  at  the  beginning  of  Appendix  C.  This 
COMDECK  is  inserted  into  a  subroutine  with  the  directive  *CALL  MAIN.  Of 
course,  information  can  also  be  passed  from  one  user  subroutine  to  another 
through  labeled  COMMON.  The  parameters  in  the  following  list  are  commonly 
needed  in  the  user  subroutines. 

NMESH  The  number  of  mesh  intervals  in  the  length. 

LAYER  The  number  of  layers  in  a  cross  section. 

NIB  =  2  The  index  of  the  mesh  point  at  end  1. 

N2B  =  NIB  +  NMESH  The  index  of  the  mesh  point  at  end  2. 

N1V  and  N2V  The  range  of  indices  of  mesh  points  that  move. 

4 . 1  INGEOM 


Subroutine  INGEOM  is  called  from  START.  It  supplies  the  initial 
geometry  and  some  control  parameters.  The  parameters  supplied  are  listed 
below. 

IB  This  is  the  dominant  control  parameter.  It  should  be 

set  to  1  or  0.  if  IB  >  0,  a  beam  is  assumed.  If  IB  <  0, 
it  is  not. 
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RADIUS  This  parameter  separates  radial  symmetry  from  slab  symmetry. 

If  RADIUS  >  0,  radial  symmetry  is  assumed  and  the  parameter 
IS  is  set  to  1  in  START.  Otherwise,  IS  is  set  to  2  in  START. 
(RADIUS  is  superseded  by  IB.  If  IB  >  0,  RADIUS  is  set  to  0 
in  START.) 

DETA1  (A?1)  This  parameter  is  a  multiplier  for  energy  computations. 

It  should  be  the  width  of  the  reference  surface.  (If  IB  =  1, 
DETA1  is  set  to  1.0  in  START.  This  makes  the  concept  of 
surface  area  and  beam  length  interchangeable.) 

DETA2  (AC)  This  is  the  distance  between  mesh  points  in  the  material 

coordinate  C  (C2  of  Ref.  1,  and  n2  of  Ref.  2).  (The  stability 
criteria  assume  that  DETA2  is  the  minimum  distance  between 
mesh  points,  and  for  some  other  coding  it  is  tacitly  assumed 
that  C  is  arc  length  from  end  1.) 

R(N) ,  Z(N)  (NIB  <  N  <  N2B)  The  position  vector- of  the  N'th  mesh 
point  on  the  reference  curve  is  rN  =  R(N)  i^  +  Z(N) 

The  point  r^  corresponds  to  the  material  coordinate  point 

C  =  (N  -  1)  DETA2. 

ETAD2  and  ETAG2(I),  I=1,NSTRN  The  previously  read  parameters  ETAD2 
and  ETAG2(I)  may  need  to  be  converted  to  the  same  units  as 
DETA2 . 

I GAUSS  Set  IGAUSS  =  1  for  Gaussian  integration  (see  Section  3.4). 

Extra  Parameters  for  Beams 

If  IB  =  1,  the  cross  section  of  the  beam  must  be  described.  It  is 
defined  as  follows: 

ZU  (;  )  The  distance  from  the  reference  surface  to  the  upper 
surface  of  the  beam. 

ZETA(K),  K=l, LAYER  (^)  These  are  the  positions  of  integration 

stations  through  the  thickness  relative  to  the  reference 
surface,  the  neutral  axis  for  bending  of  the  beam. 

(ZETA(l)  <  0,  ZETA( LAYER)  >  0) 

W(K),  K=l, LAYER  (W^)  The  cross-sectional  area  associated  with 
ZETA(K) . 

It  is  assumed  that 

W.  =  Area  of  the  beam's  cross  section, 
k  K 

£  ^k^2  Wk  =  Moment  inertia. 


E  Wk  "  0 


(i.e.  x,  =  zero  on  the  neutral  surface). 


4.1.1  INGEOM  for  a  Flat  Plate  or  Beam.  This  is  the  INGEOM  stored  with 
RPSL1D  (see  Appendix  C) . 

If  the  run  is  for  slab  symmetric  motion  of  an  initially  flat  plate, 
there  is  only  one' input  card.  The  input  is  simply  length,  SLABL,  and  one- 
half  the  width,  SLABW,  as  it  is  for  a  similar  standard  REPSIL  subroutine. 
First,  the  subroutine  reads  SLABL,  SLABW,  and  IB  with  FORMAT  (2E10.4,I5). 

The  program  sets 

RADIUS  =  0.0  , 

DETA1  =  2.0  *  SLABW  , 

DETA2  =  SLABL/FLOAT(NMESH)  , 

R(N)  =  0.0  / 

>  NIB  <  N  <  N2B  , 

Z(N)  =  (N  -  NIB)  DETA2  ) 

IGAUSS  =  1. 

If  IB  <  0,  the  responding  surface  is  an  initially  flat  plate.  The  sub¬ 
routine  is  finished. 

If  IB  >  0,  data  describing  the  cross  section  of  the  beam  is  needed. 
First,  the  subroutine  reads  ZU  with  FORMAT  (3E10.2)  and  prints  SLABL  and  ZU. 
It  then  reads  WIDTOK,  DZETAK,  and  ZETA(K)  with  FORMAT  (3E10.2),  sets  W(K)  = 
WIDTHK  *  DZETAK,  and  writes  K,  WIDTHK,  DZETAK,  ZETA(K) ,  and  W(K)  for 
K  =  1, 2, . . ., LAYER. 

The  following  table  is  an  example  for  a  6-layer  subdivision  of  the 
cross  section  of  a  solid,  right-circular,  cylindrical  rod  of  unit  radius: 


K 

WIDTHK 

DZETAK 

ZETA(K) 

1 

1.032494026 

1/3 

-0.8253794582 

2 

1.717575183 

1/3 

-0.4952276749 

3 

1.962319769 

1/3 

-0.1650758916 

4 

1.962319769 

1/3 

0.1650758916 

5 

1.717575183 

1/3 

0.4952276749 

6 

1.032494026 

1/3 

0.8253794582 

To  produce  this  table  we  first  arbitrarily  chose  layers  of  equal  thickness, 
DZETAK  =  1/3.  Then,  we  found  WIDTHK  so  that  W(K)  =  DZETAK  *  WIDTHK  would 
equal  the  area  of  the  circular  cross  section  in  that  layer.  Then  we  set 
Z „  =  K/3  -  7/6  as  a  tentative  location  in  the  center  of  each  layer.  Finally 
we  found  ZETA(K)  =  c  ZR  where  c2  ^  (Z^)  2  W(K)  =  tt/4,  the  area  moment  of 


inertia  for  a  circle  of  unit  radius. 


Similarly,  for  a  rectangular  beam  with  unit  thickness  and  width  w: 


1 


u 


L 


« 


t 


K 

WIDTHK 

DZETAK 

ZETA(K) 

ZK 

1 

w 

1/6 

-0.4225771 

-5/12 

2 

w 

1/6 

-0.2535462 

-3/12 

3 

w 

1/6 

-0.0845154 

-1/12 

4 

w 

1/6 

0.0845154 

1/12 

5 

w 

1/6 

0.2535462 

3/12 

6 

w 

1/6 

0.4225771 

5/12 

(Choosing  integration  stations  and  weights  for  Gaussian  integration,  see 
Section  3.4,  would  probably  be  better  for  a  rectangular  beam.) 

The  choice  of  layers  of  equal  thickness  is  convenient  and  reasonable 
for  these  two  examples,  but  not  for  all  cases.  With  an  I-beam,  for  example, 
it  would  be  more  reasonable,  and  convenient,  to  approximate  each  flange  with 
a  layer  and  put  four  layers  in  the  web.  If  the  beam  is  not  symmetric  ver¬ 
tically  it  may  be  difficult  to  assign  values  that-  produce  the  correct  area 
moment  of  inertia  without  altering  the  area,  or  the  location  of  the  centroid. 

4.1.2  INGEOM  for  a  Cylindrical  Shell.  The  subroutine  INGEOM  for  a 
cylinder  is  very  simple  (see  the  tabulation  in  Appendix  F) .  The  input  is 
CYLL,  RADIUS  with  FORMAT  (2E12.6),  where  CYLL  is  the  cylinder  length  and 
RADIUS  is  the  radius  to  the  midsurface.  The  subroutine  sets 

DETA1  =  2  tt  RADIUS 
DETA2  =  CYLL/ FLOAT (NMESH) 

IB  =  0 

R(N)  =  RADIUS 

Z(N)  =  (N  -  2)  *  DETA2 

If  Gaussian  integration  is  desired,  IGAUSS  =  1. 

4.2  PRFSS 

Subroutine  PRESS  supplies  the  pressure,  P(N),  at  each  mesh  point  that 
moves  (N1V  <  N  <  N2V) .  (this  is  actually  the  pressure  difference  for  shells 
and  the  normal  force  per  unit  length  for  beams,  but  we  will  continue  to  call 
it  pressure.)  The  sign  of  P(N)  was  chosen  so  that  a  positive  P(N)  tends  to 
crush  a  cylinder  or  push  a  flat  plate  down. 

Most  of  the  PRESS  subroutines  for  RPSL1D  have  been  in  two  sections. 

The  first  section  reads  input  and  makes  other  preliminary  calculations  the 
first  time  the  subroutine  is  called.  The  other  part  supplies  the  P(N)  each 
time  the  subroutine  is  entered. 

PRESS  is  called  in  the  initial  portion  of  RPSL1D  if  LOAD  =0.  It  is 
called  at  the  start  of  each  succeeding  time  cycle  in  RPSL1D  if  LPRESS  >  NCYCLE. 
Also,  as  long  as  LPRESS  >  NCYCLE,  the  P(N)  are  multiplied  by  a^  in  subroutine 
DGEOM.  This  multiplication  changes  the  P(N)  from  pressure  to  force  per  unit 
initial  area.  If  the  program  is  continued  beyond  cycle  NCYCLE  =  LPRESS,  the 
values  of  P(N)  do  not  change.  This  means  the  force  per  unit  initial  area  is 
fixed,  the  pressure  varies  with  the  change  in  area.  (At  the  initiation  of 
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damping,  when  NCYCLE  =  MDAMP*  subroutine  DAMP  sets  LPRESS  and  all  P(N)  to 
zero.  This  must  be  changed  if  the  program  is  to  continue  with  both  pressure 
and  damping.)  NCYCLE  is  the  count  of  time  cycles.  LOAD,  LPRESS,  and  MDAMP 
are  input  parameters.  If  LOAD  =  0  and  LPRESS  =  0,  PRESS  will  never  be  called 
and  all  P(N)  will  be  zero. 

One  of  the  weak  points  of  using  deformation  codes  to  simulate  reality 
is  in  accurately  describing  dynamic  pressure  loading.  It  is  difficult  to 
obtain  accurate,  repeatable  measurements  for  enough  points  even  with  well- 
controlled  experiments.  Such  measurements  tend  to  be  very  oscillatory  due 
to  both  reality  and  instrumentation.  Even  assuming  a  completely  accurate 
pressure  record  at  some  point,  the  effective  pressure  on  a  thin  shell 
would  be  different  because  of  its  rapid  reaction.  The  best  we  can  hope  for 
is  an  approximate  simulation  that  will  lead  to  useful  results.  Three  ex¬ 
amples  are  discussed. 

4.2.1  PRESS  for  Constant  Pressure.  This  is  the  PRESS  subroutine  with 
the  program  in  Appendix  C.  In  the  initial  entrance  for  each  run,  the  con- 
stand  pressure,  P0,  is  read  with  FORMAT  (E12.6).  During  each  entrance  it 
sets  P(N)  =  PO  for  N1V  <  N  <  N2V. 

4.2.2  PRESS  for  Pressure  as  a  Function  of  Time.  This  subroutine, 
called  LINPRS,  is  tabulated  with  Example  1  in  Appendix  E.  Pressure  is  a 
tabular  function  of  time.  On  the  first  entrance,  this  subroutine  reads 
pairs  of  time  and  pressure  data  into  TPR(I)  and  PPR(I)  with  FORMAT  (2E10.3) 
until  TPR(I)  <  0.  Then  IPRESS,  the  number  of  data  points  in  the  table,  is 
set  to  1-1.  (It  is  assumed  that  IPRESS  <  50,  TPR(l)  =  0.0,  TPR(I+1)  >  TPR(I), 
and  TPR (IPRESS)  is  larger  than  any  time  to  be  reached. 

Linear  interpolation  is  used  at  every  entry  to  compute  the  pressui.  PO, 
corresponding  to  time,  TIME.  The  program  then  sets  P(N)  -  F0  for  N1V  <  N  <  N2V. 

4.2.3  PRESS  for  Two  Phase  Pressure  Decay  in  a  Cyj  IrJer.  This  sub¬ 
routine,  labeled  PRESS  of  3/10/76,  is  listed  in  Appendix  F.  It  simulates 
the  pressure  in  a  cylinder  from  an  explosion  at  its  center.  This  rather 
involved  subroutine  is  included  because  it  demonstrates  a  number  of  the  com¬ 
mon  features  of  pressure  subroutines  that  simulate  the  loading  from  explo¬ 
sions.  The  unbalanced  pressure  P^  at  point  N  is  zero  until  the  shock  front 

arrives  at  time  TA^  with  peak  pressure  P0^.  The  pressure  then  decays  rapidly. 

In  this  subroutine  it  decays  rapidly,  exponentially  until  it  reaches  a  pres¬ 
sure  P  which  we  call  the  quasi-static  pressure.  It  then  continues  to  decay 
exponentially  at  a  much  slower  rate. 

PN  =  0.0  if  t  <  TAn  , 

PN  =  P0N  e'a(t  "  TV  if  TAjj  <  t  <  TBn  ,  (70) 

PN  =  P  e"Rft  “  TV  if  TBN  <  t  , 

where  a  and  6  are  the  decay  factors  and  P  is  the  quasi-static  pressure  which 

is  reached  at  point  N  at  time  TBN  found  from 


P  =  P0M  e'aCTBN  "  TV  . 

N 


An  outline  of  the  subroutine  is  given  below.  First,  the  following 
input  cards  are  read  on  the  initial  entry. 


Content 


FORMAT 


P,  6,  (PID(I) ,  1=1,6) 

CZ,  CT,  CP,  TD 
a,  — ISM 

(ZSI(I),  TSI(I),  PSI(I),  1=1, ISM) 


(2E10.3,6A10) 

(4E10.3) 

(E10.3, 10X, 110) 
(3E10.3) 


Here, 


P  is  the  quasi- static  pressure. 

$  is  the  quasi-static  decay  factor. 

PID(I)  is  an  alphanumeric  title. 

CZ,  CT,  CP,  and  TD  are  conversion  factors  to  transform  data  to  the 
desired  units  (see  below) . 

a  is  the  shock  pressure  decay  factor. 

ISM  is  the  number  of  entries  in  the  following  table  (ISM  <  50) . 
ZSI(I)  is  the  distance  from  the  middle,  end  1,  of  the  cylinder. 
TSI(I)  is  the  shock  arrival  time  at  ZSI(I). 

PSI(I)  is  the  peak  shock  pressure  at  ZSI(I). 

These  data  cards  are  read,  and  their  images  printed,  during  the  initial 
entrance.  The  tabular  data  is  then  converted  to  the  desired  units  with 
the  replacement  formulas 


ZSI(I)  =  CZ 

*  ZSI(I)  , 

(71) 

TSI(I)  =  CT 

*  (TSI(I)  -  TD)  , 

(72) 

PSI(I)  =  CP 

*  PSI(I)  . 

(73) 

This  converted  table  is  also  printed.  For  N1V  <  N  <  N2V,  TS(N)  =  TA^^  and 

PS(N)  =  P0N  are  found  by  linear  interpolation  at  distance  Z(N)  =  (N  -  2)  *  DETA2. 

(Here  is  an  example  where  the  material  parameter  £  is  assumed  to  be  length 
from  end  1.  Symmetry  is  assumed  at  end  1  where  N  =  NIB  =  N1V  =2.) 

TSB(N)  =  TB^  is  found  from 


TSB(N)  =  TS(N)  +  {ln(PS(N)/P)}/ct  . 


(74) 


The  program  then  forms  PSB(N)  and  transforms  PS(N)  to  a  more  convenient 
form,  as  follows: 
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PS(N)  =  PS(N)  e  U 

PSB(N)  =  P  e6756^  .  (7 

This  completes  the  preparation.  For  every  entry,  the  subroutine  com¬ 


putes  e  at  and  e" ^  and  then,  for  N1V  <  N  <  N2V, 
P(N)  =  0.0 


if  t  <  TS(N)  , 


P(N)  =  -e"at  PS (N)  =  -P0X,  e  “a(t  "  T/V  if  TS(N)  <  t  <  TSB(N)  ,(77) 

N 

P(N)  =  -e‘et  PSB(N)  =  -P  e'^1  _  TV  if  TSB(N)  <  t  . 

4.3  INVEL 

Subroutine  INVEL  supplies  the  initial  velocity  of  moving  mesh  points. 
The  velocity  of  all  mesh  points  is  preset  to  zero  before  INVEL  is  called. 

This  subroutine  sets 

DR(N)  =  R  ,  (78) 

DZ(N)  =  Z  ,  (79) 

where  the  velocity  at  point  N  is 

v  =  ft  i_  +  i  i,  .  (8°) 

INVEL  is  called  during  the  initiation  sequence  from  the  main  program,  RPSL1D, 
if  parameter  LOAD  is  less  than  or  equal  to  zero.  The  components  of  initial 
velocity  are  changed  to  displacement  increments  by  multiplication  with  At  in 
RPSL1D. 

4.3.1  INVEL  for  Normal  Velocity  at  Specified  Points.  The  INVEL  sub¬ 
routine  catalogued  with  RPSL1D  (see  Appendix  C)  is  a  one-dimensional  version 
of  the  INVEL  tabulated  in  the  User's  Manual  (Reference  2).  This  subroutine 
assumes  all  initial  velocities  are  normal  to  the  surface.  The  input  cards 


Cards 

,NI,NF,VR,NV 
~,  N,  V 


FORMAT 

(10X, 2I5,E12.6,I5) 
(5X,I5,E12.6) 


NI  and  NF  are  the  minimum  and  maximum  mesh  numbers  for  an  array  of  points  to 
be  given  initial  velocity  VR.  NV  is  the  number  of  mesh  points  to  receive 
initial  velocity  different  from  zero  or  VR.  N  is  the  mesh  number  for  a  point 
receiving  velocity  V. 

Note!  This  subroutine  was  written  to  use  input  cards  for  the 
standard  REPSIL.  The  mesh  numbers  NI,  NF,  and  N  are  all  relative  to 
the  mesh  point  at  end  1  being  number  1.  In  RPSL1D,  the  mesh  point  at 


end  1  is  numbered  NIB.  (Subroutine  START  sets  NIB  =  2.)  Therefore,  this 
INVEL  subroutine  adds  NIB-1  (i.e.  one)  to  all  input  mesh  numbers. 


The  first  input  card  is  read  and  the  program  sets 

DR(N)  =  -SNR(N)  *  VR  ,  (81) 

DZ(N)  =  -SNK(N)  *  VR  ,  (82) 

for  NI  +  NIB-1  <  N  <  NF  +  NIB-1,  where  the  unit  normal  at  mesh  point  N  is 

n  =  SNR(N)  i3  +  SNK(N)  i2  . 

(Notice  the  reversal  in  sign.  VR  is  assumed  to  be  an  inwardly  directed 
velocity,  possibly  from  impulsive  loading,  on  the  outside  surface.) 

The  program  then  reads  NV  cards  of  the  second  type  (NV  may  be  zero) , 
replaces  N  by  N  +  NIB-1,  and  sets 


DR(N)  =  -SNR(N)  *  V  ,  (83) 

DZ(N)  =  -SNK(N)  *  V  .  (84) 

The  subroutine  prints  the  initial  velocities  assigned.  These  initial 
velocities  are  later  converted  to  displacement  increments  by  multiplying 
them  by  At. 


4.3.2  INVEL  for  Constant  Lateral  Velocity.  This  subroutine  has  been 
used  to  simulate  a  rod  striking  a  wall.  The  changes  in  the  catalogued 
INVEL  to  produce  this  are  tabulated  with  Example  2  in  Appendix  E. 

The  constant  velocity,  VR,  is  read  with  effective  FORMAT  (20X,E12.6) 
and  the  subroutine  sets 


DZ(N)  =  -VR 


(85) 


for  N1V  <  N  <  N2V. 

4.3.3  INVEL  for  Beam  Vibration  in  the  Fundamental  Transverse  Mode. 

A  number  of  INVEL  subroutines  have  been  coded  for  beams.  They  were  either 
very  simple  or  very  specialized.  As  examples  of  simple  codes  for  initially 
flat  beams,  consider  the  following  for  the  fundamental  mode  of  free  trans¬ 
verse  vibrations.  The  transverse  velocity  of  any  point  x  at  time  t  may  be 
written 


R(x,t)  =  A  X(x)  cos(wt) 
for  these  examples.  The  velocity  at  time  t  =  0  is 

R(x,0)  =  A  X(x)  . 


(86) 


(87) 
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Here  co  is  the  fundamental  frequency,  X(x)  is  a  function  of  the  end 
conditions,  and  the  parameter  A  is  arbitrary.  We  assume  evenly  spaced 
points  to  find  x.  (These  formulas  are  similar  to  those  given  by  Enrico 
Volerra,  E.  C.  Zachmanoglou,  Dynamics  of  Vibration,  Charles  E.  Merrill 


Books,  Inc.,  Columbus,  Ohio.)  ~  " 

Cantilever  Beam: 

x  =  1.875104069  (N  -  N1B)/(N2B  -  NIB)  ,  (88) 

DR(N)  =  A (cosh(x)  -  cos(x)  -  0.7340955137  [  sinh(x)  -  sin(x)] }  .  (89) 

Simply  Supported  Beam: 

x  =  n(N  -  NIB) / (N2B  -  NIB)  ,  (90) 

DR(N)  =  A  sin(x)  .  (91) 

Clamped  Beam: 

x  =  4.730040744862  (N  -  N1B)/(N2B  -  NIB)  ,  (92) 

DR(N)  =  A  (cosh(x)  -  cos(x)  -  0.9825022145  [  sinh(x)  -  sin(x)] }  .  (93) 

4.4  INNORM 


Subroutine  INNORM  is  a  new  user-supplied  subroutine  that  assigns  the 
unit  normal  at  clamped  ends.  It  is  called  in  START  immediately  after  the 
call  of  INGEOM.  The  version  of  INNORM  catalogued  with  the  program  in 
Appendix  C  is  the  only  one  used  up  to  now.  It  produces  the  correct  normal 
if  the  shape  of  the  reference  curve  at  the  end  can  be  written  as  a  quadratic 
equation  in  E,  through  the  last  three  end  points. 

The  components  of  the  unit  normal,  ji,  are  defined  in  terms  of  partial 
derivatives  of  the  position  vector  r. 


£  =  R  i3  +  z  i2  , 

R2  =  3R/3£  , 

Z2  =  3Z/35  , 

D  =  { (R2) 2  +  (Z2)2}**  , 
SNR(N)  =  Z2/D  , 

SNK(N)  =  -R2/D  , 
n  =  SNR(N)  i  +  SNK(N)  i2  . 


(94) 


(95) 

(96) 


(97) 

(98) 

(99) 
(100) 


The  partial  derivatives,  R2  and  Z2,  at  the  clamped  ends  are  approximated  by 
non-central  differences  through  three  end  points.  N  is  NIB  or  N2B,  depending 
on  which  end  is  clamped. 

3R I 


(e.g. 


<\  r 


=  {3R(N2B)  -  4R(N2B  -  1)  +  R(N2B  -  2)}/(2AO  .) 


(101 


4.5  ENDFRC 


r 


Subroutine  ENDFRC  is  a  user  subroutine  that  supplies  the  forces  and 
moments  imposed  at  the  ends  of  a  beam.  It  is  really  a  part  of  the  option 
called  APLFRC  (Section  5.4)  which  was  programmed  to  utilize  this  type  of 
loading.  The  subroutine  is  not  needed  unless  the  APLFRC  option  is  included 
in  the  program  and  either  IBCE1  =  4  or  IBCE2  =  4.  With  the  APLFRC  option, 
IBCE1  =  4  signals  an  applied  force  and  moment  at  end  1  and  IBCE2  =  4  signals 
an  applied  force  and  moment  at  end  2.  A  free  end  is  a  special  case  for 
which  the  applied  force  and  moment  are  zero. 

Subroutine  ENDFRC  is  called  from  subroutine  RESULT  when  N  =  NIB  if 
IBCE1  =  4,  or  when  N  =  N2B  if  INCE2  =  4  and  IBCE1  ¥=  4.  Since  RESULT  is  not 
called  in  the  initial  sequence,  there  are  no  initial  forces  on  the  ends. 

When  ENDFRC  is  called,  the  force  components  and  moments  are  supplied. 
The  force  at  end  1  is 


EFR1  i,  +  EFZ1  i- 
~3  ~2 


(102) 


The  moment  is  EMI.  Similarly  at  end  2,  the  force  components  and  moment  are 
EFR2,  EFZ2,  and  EM2,  respectively.  Figure  4.5  and  the  following  list  both 
explain  the  sign  conventions. 


EFR1 

> 

0 

decreases 

R(N1B)  , 

EFZ1 

> 

0 

decreases 

Z(N1B)  , 

EM  1 

> 

0 

decreases 

R(N1B)  and 

EFR2 

> 

0 

increases 

R(N2B)  , 

EFZ2 

> 

0 

increases 

Z(N2B)  , 

EM  2 

> 

0 

increases 

R(N2B  -  1) 

EFZ 


EFZ2 


Figure  4.5  Sign  Convention  for  Applied  Forces  on  the  End  of  a  Beam. 
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The  program  outside  ENDFRC  does  not  change  these  terms.  Constant  values 
may  be  inserted  in  the  initial  entry  and  need  not  be  reset.  It  does  no  harm 
to  insert  values  at  an  end  that  is  not  forced  (IBCE1  #  4  or  IBCE2  #  4);  the 
program  ignores  them. 

4.5.1  ENDFRC  for  One  or  Both  Ends  Free.  The  initial  entry  into  this 
ENDFRC  sets  EFR1,  EFZ1,  EMI,  EFR2,  EFZ2,  and  EM2  to  zero.  These  are  never 
changed.  The  APLFRC  option  then  imposes  a  free  end  at  end  1  if  IBCE1  =  4, 
and/or  a  free  end  at  end  2  if  IBCE  =  4.  This  subroutine  is  tabulated  in 
Appendix  F. 

4.5.2  ENDFRC  for  Tabular  Forces.  This  ENDFRC  subroutine  assumes  that 
end  2  is  free  and  the  forces  on  end  1  are  supplied  by  tables.  On  the  initial 
entry  this  subroutine  reads  three  tables.  Each  has  a  heading  giving  the 


number  of  entries  and  an  alphanumeric  title. 

Cards  FORMAT 

NFZPTS,  (LABELV(I) ,  1=1,7)  (I10,7A10) 

(EFZT(I),  TMFZ(I) ,  1=1, NFZPTS)  (2E15.7) 

NFTHS,  (LABELV(I),  1=1,7)  (I10,7A10) 

(EFTHT(I) ,  TMTH(I) ,  1=1, NFTHS)  (2E15.7) 

NFMUS,  (LABELV(I) ,  1=1,7)  (I10.7A10) 

(EFMUT(I),  TMMU(I) ,  1=1, NFMUS)  (2E15.7) 


The  first  table  is  the  horizontal  force  component,  EFZ1,  as  a  function  of 
time.  The  other  two  tables  are  of  parameters  6  and  y,  respectively,  as 
functions  of  time. 

The  initial  entrance  into  the  subroutine  reads  and  prints  the  input  and 
sets  EFZ2  =  EFR2  =  AM2  =  0.  On  each  entry,  the  subroutine  finds  EFZ1,  6, 
and  y  as  functions  of  time  by  linear  interpolation.  It  then  sets 

EFR1  =  yEFZl  ,  (103) 

EMI  =  -ZU  (sin(e)  EFZ1  +  [1  -  cos(6)]  EFR1 }  .  (104) 

4.5.3  ENDFRC  for  Tabular  Forces  (Modified) .  This  subroutine  is  like 
the  previous  one  except  that  the  force  on  end  1  is  assumed  relative  to  the 
tangent,  ;r,  at  end  1, 


X  =  SNK(NIB)  i3  -  SNR(NIB)  i2  ,  (105) 

rather  than  the  horizontal  direction,  ij.  The  force  components  for  that 
subroutine  are  replaced  by 


EFR1  =  EFR1  *  SNR(NIB)  -  EFZ1  *  SNK(NIB)  , 
EFZ1  =  EFR1  *  SNK(NIB)  +  EFZ1  *  SNR (NIB)  . 


(106) 

U07) 


A  slightly  modified  version  of  this  is  tabulated  with  example  2  in 
Appendix  E. 


5.  OPTIONAL  CHANGES 

As  has  been  stated  before,  it  seems  that  any  particular  project  for 
RPSL1D  requires  some  changes  in  the  program.  Ideally,  these  changes  are 
restricted  to  the  user  subroutines,  discussed  in  the  previous  section, 
which  describe  the  geometry  and  loading.  In  this  section  we  will  briefly 
discuss  the  seven  sets  of  optional  coding  that  are  used  with  the  two  examples 
in  Appendix  E.  They  are  of  varying  utility.  The  plotting  option  PLOTP  is 
usually  wanted.  It  was  not  added  to  the  catalogued  program  because  it  is 
very  short  and  usually  requires  some  additional  coding  to  select  the  functions 
of  time  to  be  plotted.  The  APLFRC  option  which  permits  loading  the  ends  of 
beams  by  force  and  moment  has  been  used  a  number  of  times.  If  it  is  used, 
the  option  SHR3/1  must  be  included  and  the  energy  computation  EAPFRC  should 
be  added. 

It  is  hoped  that  these  options  will  be  useful  and  serve  as  a  guide  for 
introducing  others  in  the  future. 

5.1  PLOTP  (Plot  Functions  of  Time) 

This  coding  was  originally  inserted  to  plot  the  force/unit  initial  aree 
(force/unit  initial  length,  for  beams)  at  prescribed  mesh  points  as  functions 
of  time.  However,  any  function  of  time  may  be  plotted  by  inserting  coding 
that  stores  that  function's  value  into  P(N)  (N2B  <  N  <  104)  each  cycle,  and 
requesting  a  plot  for  that  P(N)  on  the  new  input  card,  card  11a.  If  the 
PLOTP  option  is  in  the  deck,  this  input  card  must  be  included  (see  Section  6) . 
It  selects  NNPE  (0  <  NNPE  <  9)  values  of  N  for  which  plots  of  P(N)  are  to  be 
made. 

The  changes  to  RPSL1D  are  rather  trivial:  some  new  terms  are  inserted 
into  COMMON  and  MAIN,  the  new  input  card  is  read  and  its  image  is  printed 
in  START,  and  the  requested  P(N)  are  stored  in  the  plot  data  file  in  PDATA. 
More  extensive  changes  were  required  in  the  plotting  program.  These  changes 
are  included  in  the  plotting  program  stored  in  file  RPSL1DPL0T  except  for 
two  orders  (the  optional  coding)  which  read  from  the  plot  data  file  and  acti¬ 
vate  P(N)  plotting.  This  option  is  used  with  both  examples  in  Appendix  E. 

(In  example  1  the  total  large  core  memory  requested  for  plotting  is  greater 
than  131,071  words.  Hence,  the  LCM  =  I  parameter  is  specified  on  the  FTN 
control  card.) 

5.2  SHR3/1  (Compute  Moment,  Shear  and  Axial  Forces 

This  option  is  for  beams  only.  It  computes  the  moment,  M,  at  each 
mesh  point  and  computes  the  shear  force,  V,  and  the  axial  force,  Q,  between 
mesh  points.  The  following  equations  were  extracted  from  Reference  1  by 
Dr.  Santiago. 


(108) 

(109) 


1 


5 


h 


ts 


Is 


-22 

M  =  a  M 

n  „  /A22  .2  £22- 

Q  =  a  (Q  -  b2  M  ) 


V  =  (l/a^)  3(a  M22/9C2)  (110) 

The  RPSL1D  analog  of  these  equations  is 

AMS(N)  =  A22  *  F22  (HI) 

QS(N)  =  A22  *  (Q22  -  BM22  *  F22)  (H2) 

VS(N)  =  (AMS(N+1)  -  AMS(N)  }/(DETA2  *  SRA)  .  (113) 


For  convenience,  the  three  new  arrays  are  put  in  COMMON  in  COMDECK  MAIN. 
The  rest  of  this  option  is  in  subroutine  RESULT.  A  table  of  N,  AMS(N),  VS(N), 
QS(N),  A22,  B22,  Q22,  and  F22  may  be  printed  at  cycles  prescribed  by  coding 
inserted  with  the  option. 

5.3  MSQSVS  (Plot  Moment,  Shear  and  Axial  Forces) 

This  option  was  inserted  to  plot  the  moment,  shear  force,  and  axial 
force  computed  with  option  SHR3/1  on  the  same  cycles  that  cross-sectional 
plots  are  made  (i.e.  whenever  NCYCLE  =  NC3DP(I)) .  The  only  change  in  RPSL1D 
is  one  statement  in  PDATA  that  stores  the  three  arrays  AMS(N),  QS(N),  and 
VS(N)  on  the  plot  data  file.  The  bulk  of  this  option  is  inserted  into  the 
plotting  program  catalogued  in  file  RPSL1DPL0T. 

5.4  APLFRC  (Applied  Forces  at  Beam  Ends) 

This  option  was  inserted  to  permit  an  applied  vector  force  and  moment 
at  either  or  both  ends  of  a  beam.  We  will  refer  to  this  as  a  forced  end. 

A  free  end  is  a  special  case  with  zero  applied  force  and  moment.  It  is 
assumed  that  the  option  SHR3/1  is  included  in  the  program. 

The  equation  of  motion,  given  in  Reference  1,  reduced  to  that  of  a 
beam  with  no  external  forces  is 

Ar+  =  Ar  +  I  M*22,  22  *  £*2»2  j  (At) 2/ (A^)  ,  d14) 

where 

M*22  =  aV2  2  (115) 

and 

N*2  =  a4  (Q22  a2  *  T^2  M22  n)  (116) 

are  vectors  defining  the  bending  and  stress  resultant  terms.  The  notation 
X, _  and  X,-,  denotes  the  partial  derivatives  9X/3^and  92X/(95)2.  The  mass 

per  initial  unit  length  for  a  beam  is  given  by  A  =  A^p  Z-W^  . 
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From  the  equations  for  shear  force,  axial  force,  and  moment, 
respectively. 


v  =  a5*  (M22,2  +  2r222  M22  =  (a  M22),2/a%  , 

(117) 

Q  =  a  (Q22  -  b2  M22)  , 

(118) 

-22 

M  =  a  M  , 

(119) 

and  the  relations 

U  2  h 

a >2  =  r22  a  • 

(120) 

.2 

2  =  ~°2  ~2  * 

(121) 

a2  =  a  I  . 

(122) 

which  are  also  taken  from  Reference  1,  we  can  transform  the  vector 

(m*22,22  +  n*2,2) 

(123) 

into  the  form 

(v  n  +  Q  t) ,  2 

(124) 

(The  shear  force  and  axial  force,  V  and  Q,  are  computed  at  all  midmesh  points, 
and  the  moment,  M,  is  computed  at  all  mesh  points,  except  at  forced  ends,  by 
option  SHR3/1.)  Vector  n  is  the  unit  normal  to  the  reference  surface  and  x 
is  the  unit  tangent. 

n  =  n,  i_  +  n.  i_  . 

~  k  ~2  r  ~3 

(125) 

~  =  nr  ~2  "  nk  ~3  ' 

(126) 

The  force  vector 

F  =  (V  n  +  Qt) 

(127) 

-22 

can  be  computed  at  all  midmesh  points  if  the  moment  M  =  a  M 
the  mesh  points. 

is  known  at 

The  APLFRC  option  assumes  a  new  user  subroutine,  ENDFRC  (see  Section  4.5), 
which  supplies  M  and  F  at  forced  ends.  In  the  following,  we  use  the  subscript 

N  to  denote  quantities  at  the  N'th  mesh  point,  N  =  NIB  at  end  1  and  N2B  at 
end  2,  and  subscript  N+^  to  denote  midmesh  quantities.  If  end  1  is  forced, 

£nib  =  EFZ1  h  +  EFR1  is  ’  Mnib  =  EM1  • 

(128) 
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If  end  2  is  forced. 


For  NIB  <  N  <  N2B, 


£n2B  *  EFZ2  i2  *  EFR2  ij  •  ^26  *  EM2 


FN^  “  ^VN+Js  2n+Js  +  %+!$  ^N+!^  C130 

Then,  the  force  per  unit  length  at  N  =  NIB  and  NIB+1  if  end  1  is  forced 
and  at  N  =  N2B  and  N2B-1  if  end  2  is  forced,  is  found  by  the  difference 
formulas: 

2^N1B  =  ^~NlB+5s  “  £n1B^^A^2^  ’  ^ 

2^N1B+1  =  ^NlB+3/2  "  £n1B+5s^A^  *  ^ 

^~'2^N2B-1  =  ^£n2B-Jj  "  En2B-3/2^A^  »  ^ 

2^N2B  =  ^~N2B  "  £n2B-%^A^2^  *  ^ 

The  first  and  last  of  these  are  not  central  differences.  The  form 

(£»2^N  =  ^~*22,22  +  ~*2,2^  * 

where  only  central  differences  are  employed,  is  used  at  all  other  moving 
mesh  points  (see  Section  3.2). 

This  option  includes  minor  changes  in  MAIN,  START,  BOUNDR,  and  BOUNDU, 
and  more  extensive  changes  in  GRAD,  RESULT,  and  MOTION,  plus  the  user 
supplied  subroutine  ENDFRC.  Option  SHR3/1  must  also  be  used,  and  option 
EAPFRC  should  be. 

5.5  EAPFRC  (Energy  for  APLFRC) 


In  the  section  for  APLFRC  we  described  the  option  which  allows  forces 
to  be  applied  at  either  end  of  a  beam.  The  optional  coding  labeled  EAPFRC, 
which  computes  the  work  on  the  beam  by  the  applied  forces  and  moments,  is 
properly  a  part  of  the  APLFRC  option. 

The  work  by  a  directed  force,  F,  applied  at  a  point  that  moves  a  vector 
distance  Ar  is 

AWp  =  F  •  Ar  .  (13 

The  work  by  a  moment  of  magnitude  M  applied  while  the  normal  vector  changes 
by  an  angle  A9  is 

AWm  =  +M  A0  .  (13 

(The  sign  in  this  equation  depends  on  the  sign  convention  used.) 


’’I 


re 


This  work  must  be  compatible  with  work  from  pressure  loading  (see 
Reference  2,  pages  45  and  46),  so  an  analogous  method  is  used  in  RFSL1D. 

That  is,  the  external  work  from  time  t-At  to  time  t,  AWp(t-^At)  and  AW^t-^At), 

will  be  the  average  of  "work  increments"  computed  at  times  t-At  and  t.  We 
compute: 

AWp(t)  =  F(t)  •  (Ar(t-!2At)  +  Ar(t+%At)  }/2  ,  (138) 

AWp(t-JsAt)  =  { AWp (t- At)  +  AWp(t)  }/2  ,  (139) 

AWM(t)  =  -  M(t)  (Ae(t-^At)  +  AOft+^At) }/2  ,  (140) 

AWM(t-JsAt)  =  (AWM(t-At)  +  AWM(t)}/2  ,  (141) 


at  the  forced  ends. 


Except  for  some  preliminary  set  up,  all  the  computation  is  carried  out 
in  MOTION.  The  values  of  AWp(t-At),  AW^(t-At),  A0(t-4At),  and  the  components 

of  Ar(t-^At)  at  any  forced  end,  is  saved  from  the  previous  cycle.  The  com¬ 
ponents  of  Ar(t+^At)  are  already  stored  in  the  arrays  DR(N)  and  DZ(N)  by 
MOTION.  We  need  to  compute  A0(t+^At).  We  define  A0(t+hAt)  by 

sin  (A0(t+^At)}  i,  =  n(t)  X  n(t+At)  .  (142) 

The  unit  normal  at  time  t  is  available, 

n(t)  =  SNK(N)  i2  +  SNR(N)  i3  ,  (143) 


but  n(t+At)  must  be  computed.  This  is  computed  from 


a2  =  9r+/3C  , 

n(t+At)  =  (ax  X  a2)/|a1  X  a2|  , 


(144) 

(145) 

(146) 


where 

r+  =  r(t)  +  Ar(t+^At)  (147) 

and  a2  is  approximated  by  a  forward  or  backward  difference. 


5.6  EROD  (Erosion  at  a  Beam  End) 


This  option  depends  on  the  applied  force  option,  APLFRO.  being  in 
RPSL1D.  The  basis  for  our  model  is  to  assume  an  end  mesh  point,  denoted 
by  N  =  NIB,  which  moves  according  to  forces  on  the  true,  nearby,  eroding 
end.  Whenever  enough  erosion  has  occurred,  the  end  mesh  point  is  shifted 
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to  the  next  mesh  point.  (The  eroding  end  cannot  approach  the  midpoint  of 
the  first  mesh  too  closely,  or  the  difference  equations  at  that  end  would 
be  unreliable.) 

We  inserted  a  new  subroutine,  ERODE,  which  supplies  the  distance,  DE, 
from  the  mesh  point  NIB  to  the  eroding  end  and  advances  the  endpoint  indices 
when  necessary.  DE  is  zero  initially.  Whenever  DE  is  greater  than  AS/4, 
we  subtract  AC  from  DE  and  add  1  to  NIB,  N1V,  and  N1A  (these  are  all  indices 
for  the  end  point) .  The  difference  equations  for  the  shear  in  RESULT  and 
force/unit  length  on  the  end  point  in  MOTION  were  changed  to  reflect  the 
position  of  the  eroding  end.  The  ENDFRC  subroutine  was  altered  to  use  quad¬ 
ratic  interpolation  to  find  the  normal  at  the  eroding  end.  (This  last  change 
was  not  really  tested  since  the  rod  remained  straight  in  our  application.) 
Several  minor  inserts  were  made  to  permit  reasonable  plots. 

This  option  is  not  as  well  tested  as  other  features  of  RPSL1D.  This 
brief  discussion  is  included  because  the  option  is  used  with  Example  2 
(Appendix  E) .  The  motivation  for  Example  2  was  to  demonstrate  that  the  RPSL1D 
program  could  reproduce  experimental  strain  records  at  two  points  along  a 
long  rod  penetrator,  assuming  a  stress-strain  curve,  and  a  table  of  force 
on  the  end  and  position  of  the  end  as  functions  of  time,  all  derived  from 
the  Karmann-Taylor  theory3.  From  the  two  strain  records,  the  density,  and 
the  cross-sectional  area  of  the  rod,  one  can  derive  the  stress  and  force  for 
any  strain,  and  the  velocity  at  which  the  strain  is  propagated.  From  this 
velocity,  one  has  a  linear  relation  between  initial  time  and  position  of  the 
generating  force  for  each  strain.  Two  quite  different  tables  of  erosion  and 
force  as  functions  of  time  were  generated,  one  for  a  slow  erosion  rate  and 
the  other  for  a  rapid  erosion  rate.  Both  were  used  as  input  for  RPSL1D,  and 
both  produced  strain  records  in  fair  agreement  with  the  test  data. 

5.7  BSTRS  (Stress-Strain  Option  for  Beams) 

The  optional  coding  labelled  BSTRS  is  a  stress-strain  routine  for  beams, 
which  models  the  uniaxial  curve  more  smoothly  than  the  mechanical  sublayer 
model  by  using  more  points.  The  mechanical  sublayer  model  in  the  catalogued 
routine  is  generally  superior,  but  the  number  of  segments  in  the  polygonal 
stress-strain  curve  must  be  restricted. 

In  the  erosion  problem  discussed  in  Section  5.6,  we  wanted  to  compare 
the  experimental  strain- time  response  in  a  bar  with  output  from  RPSL1D.  The 
strain  response  from  RPSL1D  at  points  along  a  bar  agree  closely  with  those 
predicted  by  the  Karmann-Taylor  theory.  That  is,  the  velocity  at  which  a 
particular  strain  level  is  propagated  down  the  bar  is  proportional  to  the 
square  root  of  the  slope  of  the  stress-strain  curve  at  that  strain.  There¬ 
fore,  the  strain-time  records  from  RPSL1D  at  points  away  from  the  end  of  the 
rod,  in  response  to  an  increasing  force  at  the  end,  have  near  discontinuities 
corresponding  to  the  corners  of  the  polygonal  stress-strain  curve  where  there 
are  discontinuities  in  the  slope.  To  remove  these  discontinuities,  one  needs 


3 

A.D.  Gupta  and  J.D.  Wortman ,  "An  Eroding  Long  Rod  Penetrator  Model  for  Hard 
Target  Penetration , "  Proceedings  of  the  Third  ASCE/EMD  Specialty  Conference , 
September  17-19,  1979 ,  University  of  Texas,  Austin ,  Texas,  pp  714-717, 
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a  smoother  stress-strain  curve.  This  was  impractical  with  the  mechanical 
sublayer  model,  so  the  BSTRS  option  was  coded.  If  enough  points  are  used, 
the  uniaxial  loading  curve  will  be  smooth  enough.  This  option  is  not  pro¬ 
grammed  for  a  strain  rate  effect  nor  for  biaxial  loading  and  there  is  no 
obvious  simple  way  to  include  either  of  them. 

The  positive  portion  of  the  uniaxial  stress-strain  curve  is  entered  as 
a  table  of  points,  (e^,o^),  with  FORMAT  (2E15.7),  which  terminates  with 

eNEST  <  (NEST  <  100).  Then  e^g^  is  replaced  by  1000.0  (infinity)  and 

°NEST  reP^ace^  by  o^g^  As  before,  for  the  mechanical  sublayer  model, 

the  point  (e^.o^)  is  replaced  by  (SIGZ/E,  SIGZ),  the  elastic  limit.  Also, 

as  before,  these  are  not  engineering  stress  nor  strain.  If  (ee^»ae^)  are 

engineering  values. 


a. 

l 


e  .  (1  +  Jje  . )  , 
eiv  2  ei'  * 


a  . (1  +  a  . )  . 
ei^  ei' 


(148) 

(149) 


2  2 

These  correspond  to  the  mixed  components  ^2*  °2  the  strain  and  stress 

tensors  used  in  RPSL1D.  The  negative  portion  of  the  curve  is  defined  by 
a(-e)  =  -a(e)  . 


For  each  integration  station  at  each  mesh  point  and  at  each  midmesh, 
we  keep  track  of  the  strain,  e  =  e^,  and  the  mean  value  of  possible  elastic 

strain  variation,  e  ,  in  addition  to  o  = 


. 


Initially  we  set  e  =  e  =  a  =  0.0  . 

m 

Assume  we  have  e,  a,  and  e  from  time  t-At  (call  these  e”,  a”,  and  e  ”), 

m  *  '  m 

and  Ae,  the  change  in  e  from  time  t-At  to  time  t.  We  recognize  five  cases 

for  e  >  0.  (For  e  ~  <  0,  we  change  the  sign  on  e”,  o~,  e  ",  and  Ae,  and 
in  m  b  b  m  ’ 

proceed  as  for  em  >0.  At  the  end,  we  change  the  signs  on  e,  o,  and  em<) 

In  all  five  cases  we  set 

e  =  e  +  Ae  .  (150) 

Case  Conditions 

(1)  Ae  >  0 

G  <  em  +  el 

(2)  Ae  <  0 

e  >  c  -  e . 

m  1 


Results 


o=o  +  E  Ae 


c  =  e 

m  m 


o=o  +  E  Ae 


e  =  c 

m  m 


« 
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Case 


Conditions 


Results 


(3) 

(4) 

(5) 


Ae  >  0 

e  >  e  +  t. 
m  1 

Ae  <  0 

e  <  e  -  e. 


m 

e  >  -e. 


Ae  <  0 

e  <  e  -  e, 

m  1 

c  <  -e. 


a  =  a(e) 

e  =  e  -  e. 
m  l 

a  =  a(e  +  2e^)  -  2a^ 

e  =  e  +  e. 
m  1 


a  =  -a(-e) 

em  =  e  +  e, 
m  1 


The  functional  value  o(e)  is  found  from  the  table  (e^,  o .)  by  linear  inter¬ 
polation. 


(We  should  mention  that  the  computation  of  strain  energy  is  as  if  un¬ 
loading  to  o  =  0  can  take  place  entirely  elastically.  That  is,  the  strain 
energy  is  assumed  proportional  to  (a)2  where  a  =  E  e.  This  assumption,  also 
used  in  REPSIL  and  the  catalogued  RPSL1D,  is  only  true  if  |a|  <  20^.) 


The  UPDATE  cards  for  this  option  include  some  changes  and  additions  to 
COMMON  in  COMDECK  MAIN,  a  few  orders  to  initiate  some  new  arrays  to  zero  ir 
the  main  program,  RPSL1D,  a  change  to  read  and  store  the  new  form  for  the 
stress-strain  table  (input  cards  7)  in  subroutine  START,  and  a  complete  re¬ 
placement  of  BMSTRS,  the  stress  subroutine  for  beams. 


6 .  INPUT 

RPSL1D  originally  used  the  same  input  as  the  standard  REPSIL  (see 
Reference  2,  pp  54-71) .  This  is  still  true  if  the  user-coded  subroutines 
are  compatible.  Table  6.1  lists  the  input  for  RPSL1D.  The  first  part  of 
this  table  is  like  table  3.1  of  Reference  2.  Variables  that  are  no  longer 
used  have  a  line  drawn  through  them.  A  brief  description  for  each  variable 
is  given  after  the  table.  More  complete  descriptions  for  some  variables  may 
be  found  in  Reference  2.  These  descriptions  are  still  valid  and  will  not  be 
repeated  here.  Units  are  indicated  in  brackets:  T  for  time,  L  for  length, 
and  F  for  force.  Limits  are  those  for  the  catalogued  RPSL1D  (Appendix  C) ; 
most  of  these  can  be  easily  increased. 

The  input  for  user  subroutines  is  listed  in  the  order  the  subroutines 
will  be  called.  Examples  of  input  are  given,  but  the  descriptions  of  the 
input  variables  for  the  user  subroutines  are  ..  *^posely  sketchy.  The  user 
must  check  these  subroutines  and  their  input.  iNGEOM  is  always  called  first 
and  usually  requires  input  from  cards.  INNORM  is  called  next;  no  input 
has  been  needed  so  far.  PRESS  is  called  next  if  LOAD  *  0.  INVEL  is  called 


next  if  LOAD  <  0.  (PRESS  will  be  called  next  if  LPRESS  >  0;  so  far  no  ver¬ 
sion  of  PRESS  has  been  coded  for  RPSL1D  that  required  input  after  the  initial 
entrance.  The  initial  entrance  should  occur  with  LOAD  =£  0,  but  it  could 
occur  here.)  ENDFRC  will  be  called  next  if  option  APLFRC  is  used.  Subroutine 
ERODE,  part  of  option  EROD,  would  be  called  last. 


CARD 

VARIABLES 

FORMAT 

1 

TITLE 

8A10 

2 

MESH,  NMESH,  LAYER,  YLDFAC 

3I5,E12.6 

3 

MAXC,  NC0NT,  WRTTE-,  DELTAT 

3I5.E12.6 

4 

IBGE1,  IBCE2,  IBCE3,  IBCE4 

1615 

S 

LOAD,  LPRESS,  MDAMP,  DAMPF,  DFACT 

315, 2E12.6 

6 

E,  FNU,  SIGZ,  RH0,  THICKN,  NSFL,  ISR 

5E12.6, 215 

7 

(SSIG(J) ,  SEPS(J) ,  DSR(J) ,  PSR(J) ,  J=1,NSFL) 

4E15.7 

8 

NPRINT,  (JCHK(I) ,  1=1,3) 

1615 

9 

NUMCY,  (NCYCH(I) ,  I=1,NUMCY) 

1615 

10 

NLPRIN,  (JCYNLP(I) ,  I=1,NLPRIN) 

1615 

11 

N3D,  (NC3DP(I) ,  1=1, N3D) 

1615 

11a 

NNPE,  (NPE(I) ,  I=1,NNPE) 

1615 

12 

ETAD1,  ETAD2,  NSTRN 

2E10 .4, 15 

13 

(STA&H-t},  ETAG2(I) ,  ANGLE (I)  ,  ANDLB(I)  , 

NETAG(I) ,  1=1, NSTRN) 

4E10.4, 15 

14 

Input  for  subroutine  INGEOM 

15 

Input  for  subroutine  INNORM 

16 

Input  for  subroutine  PRESS 

17 

Input  for  subroutine  INVEL 

18 

Input  for  subroutine  ENDFRC 

19 

Input  for  subroutine  ERODE 

Table  6.1  List  of  Input  Cards 


o  IBCE4  from  card  4  is  stored  in  position  IBCE1  and  used  as  end  1 
condition  control. 

o  Omit  card  7  if  NSFL  =  0  (or  if  NSFL  =  1  and  ISR  =  0) ,  unless  the  BSTRS 
option  is  used.  Card  7  is  read  with  a  different  form  and  FORMAT 
with  option  BSTRS. 

o  Card  11a  is  included  if,  and  only  if,  the  PLOTP  option  is  used. 

The  input  variables  in  table  6.1  are  described  below.  Additional 
descriptions  are  given  in  Chapter  3  of  Reference  2. 

Card  1  TITLE  Alphanumeric  statement  for  output  identification. 


Card  2 


Not  used. 

Number  of  meshes.  (NMESH  <  102) 

Number  of  integration  stations  through  the  thickness. 
This  may  be  the  number  of  layers  as  in  REPSIL. 


IICOU 

I'lLUIT 

NMESH 

LAYER 


YLDFAC  Parameter  controlling  substep  size  for  plastic  yielding 
in  shells.  (2.0  is  suggested) 

Card  3  MAXC  Final  time  step.  (Total  number  of  time  steps  to  final 

time.) 

NCONT  Initial  time  step  for  run.  This  must  be  zero.  Restart 

capability  has  been  deleted  from  the  CDC  version. 

WRITE  Not  used. 

DELTAT  Suggested  time  step.  This  will  be  replaced  if  it  is 
larger  than  the  computed  stable  At  or  if  entered  as 
0.0.  [T] 

Card  4  IBCB1  Numbers  for  end  conditions.  IBCE1  and  IBCE3  are  unused. 

IBCE2  IBCE2  is  for  end  2  (C  a  maximum) .  IBCE4  is  for  end  1 

I-BGE3  (£  a  minimum).  IBCE4  is  stored  in  FORTRAN  position  IBCE1) . 

IBCE4  1  -  clamped  end, 

2  -  symmetry  end, 

3  -  hinged  end, 

4  -  free  end,  end  2  only  with  catalogued  RPSL1D, 

forced  end,  either  end  with  APLFRC  option. 

Card  5  LOAD  Parameter  which  controls  calls  of  INVEL  and  PRESS  in  the 

initiation  sequence.  (If  LPRESS  >  0,  LOAD  0) 

1  -  call  PRESS  (supplies  initial  pressure) 

0  -  call  INVEL  (supplies  initial  velocity) 

-1  -  call  PRESS  and  then  call  INVEL 

LPRESS  Subroutine  PRESS  is  called  each  time  cycle  to  supply 

pressures  in  P(N),  and  then  P(N)  is  converted  to  force/ 
initial  unit  surface  area,  until  NCYCLE  >  LPRESS.  The 
force/initial  unit  area  is  constant  after  NCYCLE  >  LPRESS. 
(For  beams,  unit  area  is  replaced  by  unit  length.) 

MDAMP  If  NCYCLE  >  MDAMP,  the  damping  procedures  are  carried 

out.  If  NCYCLE  =  MDAMP,  the  program  sets  LOAD  =  0, 

LPRESS  =  0,  and  all  P(N)  =  0. 

DAMPF  Viscous  damping  coefficient.  [FT/L3] 

DFACT  Parameter  controlling  termination  of  program  during 

damping.  (Suggest  .001) 

Card  6  E  Young's  modulus.  [F/L2] 

FNU  Poisson's  ratio,  v. 

SIGZ  Yield  stress,  oQ.  [F/L2]  (See  card  7.) 

RHO  Initial  mass  density  per  unit  volume,  p.  [FT2/L4,  (i.e. 

mass/L3) ] 

THICKN  Thickness  of  shell.  [L] 

NSFL  Number  of  changes  of  slope  in  the  polygonal  approximation 

of  the  uniaxial  loading  curve.  (NSFL  <  5) 

0  -  elastic  behavior,  no  plasticity. 

1  -  elastic  perfectly  plastic  response, 

>1  -  elastoplastic,  strain  hardening  response. 
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ISR  Strain  rate  sensitivity  control. 

0  -  no  strain  rate  effect, 

1  -  plasticity  is  strain  rate  dependent. 

Card  7  SSIG(J),  Stress,  a  [F/L2],  and  strain,  e  [L/L],  at  the  comers 

SEPS(J)  of  the  polynomial,  uniaxial  loading  curve.  (SSIG(l)  is 
replaced  by  SIGZ  and  SEPS(l)  by  SIGZ/E.)  (If  (es,as) 
are  'engineering'  strain  and  stress,  a  =  os(l  +  es) , 

e  =  es(l  +  es/2).  ) 

DSR(J),  Constants  for  strain  rate  behavior.  These  must  be 
PSR(J)  included  if  ISR  =  1. 

With  the  BSTRS  option,  cards  7  become: 

(EEB(I) ,  SSB(I) ,  I  =  1,2,...)  FORMAT(2E15.7) 

EEB(I) ,  Table  of  strain  and  stress  (reverse  order  from  above) 
SSB(I)  for  up  to  100  points.  Table  is  terminated  by 

EEB(NEST)  <  0.  SSB(l)  is  replaced  by  SIGZ,  and  EEB(l) 
by  SIGZ/E,  EEB(NEST)  by  1000,  and  SSB(NEST)  by 
SSB (NEST-1) . 


Card 

8 

NPRINT 

Number  of  cycles  between  surface  strain  prints  at  points 
defined  on  cards  13.  Maximum  deflection  and  extreme 
strains  are  also  printed. 

JCHK(I) 

Numbers  controlling  printing  at  cycles  on  cards  9. 

JCHK(l) ,  or  JCHK(2),  >  0  requests  prints  of  coordinates 
displacement  increments,  and  pressure  at  mesh  points. 

JCHK(3)  >  0  requests  prints  of  surface  normals,  and 
a  print  of  surface  strains  on  both  surfaces,  at  all 
mesh  points. 

Card 

9 

NUMCY 

Number  of  cycles  for  JCHK(I)  controlled  prints  (card  8) 
and  energy  balance  summary.  (<51) 

NCYCH(I) 

Cycles  for  prints. 

Card 

10 

NLPRIN 

Number  of  JCNLP(I)  to  follow.  (<51) 

JCNLP(I) 

Cycles  to  print  IMAT(N,K)  and  IMATM(N,K)  arrays. 

Card 

11 

N3D 

Number  of  cycles  at  which  displacement  plots  are  to  be 
made.  (<51) 

NC3DP(I) 

Cycles  for  plots. 

Card  11a  (Include  if  PLOTP  option  is  used.) 

NNPE  Number  of  plots  to  be  made.  (0  <  NNPE  <  9) 

NPE(I)  Index,  N,  of  P(N)  to  plot  as  a  function  of  time.  P(N) 

(N1V  <  N  <  N2V)  is  force  per  initial  unit  area  (initial 
length  for  beams) .  Other  functions  of  time  may  be 
plotted  by  storing  them  in  P(N)  (N2B  <  N  <  104)  and 
listing  N  on  this  card. 


Card  12  ETAD1  Not  used. 

ETAD2  Material  coordinate  of  location  at  which  displacement 
as  a  function  of  time  is  to  be  plotted. 

NSTRN  Number  of  locations  at  which  surface  strains  listed  on 

cards  13  are  to  be  plotted  and  printed  (see  card  8). 

(<  6) 

Card  13  ETAGl(I)  Unused. 

ETAG2(I)  Material  coordinate  for  surface  strains. 

ANGLE(I),  Angles  for  additional  strain  calculations.  (See  Figure 
ANGLB(I)  3.5,  Reference  2) 

NETAG(I)  Selects  surface  foT  strain. 

0  -  Surface  toward  positive  normal, 

1  -  Surface  away  from  positive  normal. 

[Sample  Input  for  User's  Subroutines.  See  USER  SUBROUTINE  section.] 

Card  14  [Subroutine  INGEOM] 

[INGEOM  for  semi-infinite  flat  plate  or  beam.  Catalogued  INGEOM] 
(SLABL,  SLABW,  IB  (2E10.5.I5)) 

SLABL  Length. 

SLABW  Multiplier  for  energy  computation,  A^1  =  2.0  *  SLABW. 

Not  used  for  beams. 

IB  IB  >  0  denotes  a  beam. 

If  IB  >  0,  also  read 

(ZU  (E10.2)) 

( (WIDTH (K) ,  DZETA(K) ,  ZETA(K) ,  K  =  1,..., LAYER)  (3E10.2)) 

ZU  Location  of  upper  surface. 

WIDTH (K)  Width  of  K'th  layer. 

DZETA(K)  Thickness  of  K'th  layer. 

ZETA(K)  Integration  station  for  K'th  layer. 

[INGEOM  for  a  right  circular  cylinder] 

(CYLL,  RADIUS  (2E12.6)) 

Card  15  [Subroutine  INNORM] 

(No  input  for  present  INNORM) 

Card  16  [Subroutine  PRESS] 

[Constant  pressure.  Catalogued  PRESS] 

(P0  (E12.6)) 

PO  Pressure.  Force  per  unit  length  for  beams. 


[Time  dependent  pressure,  LINPRS] 

( (TPR(I) ,  PPR(I),  1=1,2,...)  (2E10. 3) ) 

((-1.0)  (2E10.3)) 

[PRESS  with  tabular  data,  exponential  decay.  3/10/76] 

(PBAR,  BETA,  (PID(I),  1=1,6)  (2E10.3.6A10)) 

(CONCZ,  CONCT,  CONCP ,  TDIF  (4E10.3)) 

(ALPHAS,  NOUSE,  ISM  (2E10. 3, 110) ) 

( (ZSI (I) ,  TSI(I),  PSI(I),  1=1, ISM)  (3E10.3)) 


Card  17  [Subroutine  INVEL] 

[Catalogued  INVEL.  Like  INVEL  in  Reference  2.] 

(--,  NI,  NF,  VR,  NV  (10X, 2I5,E12.6,I5)) 

(-,  N,  V  (5X,I5,E12.6)) 

NI,  NF  Mesh  points  NI+1  through  NF+1  are  given  initial  normal 
velocity  -VR.  Mesh  point  at  end  1  is  numbered  NIB  =  2. 

VR  Velocity.  [L/T] 

NV  Number  of  cards  to  follow. 

N  Mesh  point  N+l  given  initial  normal  velocity  -V. 

V  Velocity.  [L/T] 

Card  18  [Subroutine  ENDFRC.  Needed  with  APLFRC  option.] 

[End  2  free.  Tabular  forces  on  end  1.] 

(NFZPTS,  (LABELV(I) ,  1=1,7)  (I10,7A10)) 

( (EFZT(I) ,  TMFZ(I) ,  1=1, NFZPTS)  (2E15.7)) 

(NFTHS,  (LABELV(I) ,  1=1,7)  (I10.7A10)) 

( (EFTHT(I) ,  TMTH(I) ,  1=1, NFTHS)  (2E15.7)) 

(NFMUS,  (LABELV(I) ,  1=1,7)  (I10.7A10)) 

( (EFMUT(I) ,  TMMU(I) ,  1=1, NFMUS)  (2E15.7)) 

Card  19  [Subroutine  ERODE.  Part  of  EROD  option.] 

(NTEPTS,  (LABEL(I) ,  1=1,7)  (I10.7A10)) 

((TEROD(I) ,  XEND(I) ,  1=1, NTEPTS)  (2E15.7)) 
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APPENDIX  A 


LIST  OF  EQUATIONS 


In  this  section,  we  list  the  principal  equations  used  in  the  standard 
REPSIL  and  the  corresponding  equations  used  for  axial  symmetry,  slab  symmetry, 
and  beams,  respectively.  The  notation  in  column  1  is  that  of  Reference  1, 
except  we  have  assigned  an  orthonormal  cartesian  vector  basis  (i^,  i.2>  i^) 

and  defined  corresponding  Cartesian  components  of  position,  Y1,  Y2,  and  Y3. 

The  notation  'same*  in  any  column  means  that  the  entry  is  the  same  as 
that  in  the  column  to  the  left. 
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APPENDIX  B 


DEFINITION  OF  FORTRAN  VARIABLES 


In  this  appendix  we  list  the  FORTRAN  variables  used  in  RPSL1D  as  listed 
in  Appendix  C,  plus  additional  variables  used  in  the  options  PLOTP,  SHR3/1, 
MSQSVS,  APLFRC,  EAPFRC,  EROD,  or  BSTRS.  Extra  columns  and  a  lot  of  symbology 
were  used  in  an  attempt  to  shorten  the  appendix  but  still  give  enough  infor¬ 
mation. 

Column  1  is  the  FORTRAN  name.  Closely  related  names  may  be  grouped. 
Arrays  are  indicated  by  enclosing  a  symbolic  integer  in  parentheses 
(e.g.  AMS(N) ) . 

Column  2  lists  the  symbol  used  in  this  text  or  one  of  the  references. 

If  the  symbol  is  in  braces,  {  },  it  was  first  used  in  this  report;  if  in 
brackets,  [  ],  it  was  first  used  in  Reference  2;  if  in  neither  braces  or 
brackets,  the  symbol  originated  in  Reference  1.  Examples: 

{M}  Symbol  for  moment  in  this  report. 

[A]  Coefficient  used  in  Reference  2. 

a^  Symbol  of  a  covariant  metric  component  first  used 
in  Reference  1 . 

Column  3  partially  describes  the  status  of  the  variable: 

C  denotes  a  name  in  COMMON. 

I  denotes  input. 

Sq  denotes  a  variable  that  is  rarely  changed  throughout  a 
u  computer  run,  as  opposed  to 

S  denoting  a  term  that  is  fixed  for  one  time  step,  or 

T  denoting  a  transient  value. 

D  denotes  a  dummy  argument  for  a  subroutine. 

Column  4  names  the  subroutine  which  stores  the  variable.  (If  other 
subroutines  change  or  store  these  variables,  they  will  usually  be  listed  in 
Column  6.) 

Column  5  gives  the  stage  of  development  when  the  variable  originated: 

U  denotes  a  variable  described  in  Appendix  C  of  the  User's 
Manual  (Reference  2) , 

R  denotes  a  newer  name  included  in  the  RPSL1D  file, 

any  name  is  the  option  that  first  used  the  variable. 
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I  MCI  (AW  I  T  MOTION  F.AI’FRC  Work  increment  [t,t+At],  from  moment  on  the 
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surface  in  and  £2  directions.  (Output) 


EROOP  S  EROOE  EROO  Erosion  from  initial  end  1  up  to  time  t-At. 
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NCYCHOl  C.I.S.  START  0  Cycles  for  JCIIK  and  energy  sinmiary  print.  (Input  9) 


Position  for  unwanted  input. 


2 

«l 

Cl 


C 

u 

fiC 

g 


I 


u: 

z 


8. 

C 


Z  « 


w  E 
O  •*-> 
*C  *- 


re  ^ 
c 

V.  P  4-» 

**  — •  a 
o  ^  c- 
—  o  c 
p.  c  *-. 


*4 

p. 


12 

a. 

Z 


8. 

s 


3 

P- 

C 


re  x>  re  e  x  §« 

u  «-»  o  to  s 

O  ♦-»  X  -*  4- 


o 

X 


Su 

c 

•*■*  O  »-' 


re 

3 

*3 


S5 

.oSd 

rt  O 


0> 
tJ 
*  © 

5  £ 


75 

a  £ 

a)  O 

ig 


z  *o 
o 


•  VI 
O  01 
>  T3 

8  2 
0) 


8. 

S 


o 

p. 


c 

»M 

8. 

*w 

O 

E  * 
X 


C 

CS  01 
Cl  — ■ 
I-  **s. 
CO  o 

u 

*-« 

<0  c 

•M  <M 


o 

o 

os 

=>  u: 


8  0  o 

o  o 

a.  cl  ac 

0C  M  CL  Vi  KU  O 


C/3 

c- 

oe 


ll 

ir.  c/> 


o 


H 

C/3 


I 


£ 


§ 

c/3 


a; 

C 


1 


u 
f-  o 
os  O 
<  CL 

H  ui 
C/3 


<  os 
H  u: 
c/3  w 


I-  o 

Si 

f-  U) 
C/3  w 


os 

c- 


C/3 

CO 

to 

c 

V. 

mi 

M4 

* 

a 

u 

Cm 

o 

u 

C/3  *— 1  »— ■ 


C/3 

U 


*-•  C/3 

H  U 


5  2 


ac 

Z 


»  :•  g  s 


I 


C/3 

£ 


< 

<N 

z 

< 

M* 

z 


> 

z 


o 

’p. 


z  <-> 
c. 

M- 

TJ 

Cj  / — • 

o  oo 

p. 

< 

b« 

w- 

to 

<D 

10 

to 

Qi 

•M 

TJ 

K 

w 

C 

re  4-» 

X  c 

s 

•M  (0 

C  'H 

•« 

r~  -*s 
O  C 

esc 

C 

Cm  M 

e 

• 

to 

•H 

£ 

V- 

•M 

5  8. 

00 

i 

til. 

s 

/-^ 

P. 

•M 

*-> 

3  E. 

= 

c 

re 

s 

C 

Lf 

X 

z 

4-1 

•M  W 

UJ 

A 

1 

o 

W  C 
•— « 

c 

♦-* 

C 

re 

(0 

z 

•H 

f? 

TJ 

e 

C  TJ 

M4  p 

a 

3 

z 

A 

re 

M» 

b  ' 

»H 

a> 

o 

o 

tu 

* 

5  ai 

• 

0/ 

o  01 

o 

X,  • 

H 

z 

■H  CO 

oc 

(m 

10 

o 

c 

u 

• 

w 

to 

Cm 

Cm 

b  *t5 

to 

U  II 

*  6 

0 

b 

0 

c 

M 

X 

«2 

o 

5 

e 

o 

'S 

m4 

o  ’ 

0> 

o 

C 

p. 

CM 

o  re 

u 

>1 

as 

>* 

•M 

«m  to 

01 

10 

Cm  to 

6 

cm 

z 

o 

re  w 

X 

p 

to  £ 

8. 

eo 

b 

■ 

to 

£ 

Z 

Z  to 

CO 

z  5! 

M 

CL 

4-> 

c» 

/■“S 

z  • 

X 

1 

O  4J 

Uh 

c 

P. 

«4H 

Cm 

Cm 

Cm 

Cm  to 

C  3 

0) 

er  o 

A  L 

Cm 

p 

c 

^  z 

"2 

3 

£ 

.5  .5 

i. 

V 

° 

° 

o 

*c 

u 

U  M 

U 

° 

U  /■*> 

Cm 

CM 

A 

u 

4-* 

01 

re 

£ 

O 

a> 

Cl 

x  £ 

4-»  *M 

o 

•M 

4->  *M 

0) 

M  ^ 

X 

■3 

° 

10 

l-H 

X  *H 

| 

10 

1 

E 

S 

e 

3  £ 

s  < 

a 

00 

e  > 

E 

b 

0i 

z  ^ 

•£ 

b. 

Z  p. 

z 

z 

os 

z 

Z 

£ 

z 

Z  w 

3  z 

z 

z 

-J  z 

Z 

C.  ' 

< 

> 

5 

o 

c. 


C/3 

w 

I 


C/3 

cJ* 


< 

£ 


H 

£ 


65 


STRAIN  II  3.141592653589793 .  (Too  long  for  CDC) 


g 

uj 

g 


IA  y-> 
C  «  N 

«  M  M 

U  v>  3 

m  a 

nooe 
I  c  m. 

IA  -M  «—/ 
CA  V. 

4)  4>  ^ 
Vi  4>  • 

E  MS 
IA  «M  IA 
DC  0) 

C  + 
a  4/  —• 

•M  W 

X  DC  (A 

re  co 

*c  *o  u 
s  c 


S.  cT" 


B. 

c 


o 

Ca 


<m 

Vi 

7 


U 

O 

*m 


Cm 

O 


o  V) 
re  • 

V<  Vi  /— > 

m  m  n 
c  o  ^ 

U  o  tA 

o  re  ms 
re  re  +  /— * 

«-»  M 
X  CA  w  O 


•7 

0) 


IA  (A 

CA  O 

4)  U 

Vi  V» 

4-*  IA 

(A 

*a 

•o 


I 

E 


4> 

tx 

o 


(A  0  «M  « 

re  >  ce  — « 

Vi  H  Vi  w 

V  3  V  w 

CO  4>  IA  w 


V 

V 
re 
7 
O' 
CO 


z 

h 

K 

o 

u. 

os 

2 

u 

i 

£ 

8 

UJ 

e 

to 


*■0 

CO 


CO 

« 

UJ 


Vi 

re 

X 


Vi 

re 

x 

(0 


S 

E 

*  0) 

CO 

•M 

IA 

re 

A 

V-i 

B 

IA 

4> 

*~i  Ms 

X 

IA 

re 

IA 

o 

•H 

M 

& 

IA  II 

UJ 

re 

•H 

X 

IA 

X 

C 

M  D 

•'■N* 

u 

M 

X 

re 

E 

o 

re 

0> 

re  »  *m»  m 

•M 

M 

re 

Vi 

o 

u 

v» 

IA 

IA  U 

o 

E 

CA 

•M 

M 

£X 

IA 

E  O 

Ml 

re 

re 

M 

CA 

M 

E. 

E 

X 

•H  nm  * 

CO 

E 

*7 

re 

re 

o 

M 

re  m 

X 

*x 

M 

•a 

(A 

o 

Vi  <M  E 

•a 

re 

re 

IA 

M 

O 

re 

*~3 

M  l-l  «M 

II 

•M 

•H 

re 

E. 

re 

X 

IA  w  O 

cm 

X 

cm 

hH 

E 

V 

cm 

Da 

Ms 

o 

7 

0 

X 

re  m  re 

c  re  7 

x  m  o* 

O  CO  CO 


re 

4-» 

CO 


E 

•H 

a 

•s~ 

4)  C 
6  •-* 
-J 

*  co 

K& 

V  * 

£S2 

«A  H 

-a 

O  CD 


4)  ^ 

Is 

e  x 
o  re 

O  ri 

x 

^  j 
B  CA 

re 

•H  « 

Vi  * 
re 

>  e 
re  o 


e  re 
o  *■» 

CJ  (A 


z 

X 

(A 

re 

i 

2S 

EG  g 

-  CO  K 

<A  Cm  V* 

«A  DC  CO 


w  IA 

IA  (U 

M  Vi 

B  •  *> 

4>  *->  W 

c 

O  Vi  CM 

7.  «  O 

B  X 

ore  (A 

O  M  «■» 

X  E 

4-»  7  4> 

C  IA  E 

re  o 

•f-t  -  Cu 

E 

re  o 

>  E  O 

re  o 

Vi  -H  "O 

mm  O 

Ere  x 

O  M.  •«* 

U  IA  Z 


h 

H 

h 

CO 

CO 

cc 

a: 

1 

CO 

UJ 

UJ 

& 

1 

E 

CO 

E2 

g 


£■  a 


V  . 

c 

a 

§ 


*7 

4) 

X 


§ 

§• 

o 


-o 

re 

X 


2  i 


re 

a 


re 

E 

a 

§ 


■o 

0) 

X 


E 

•M 

a 

x 

(A 

4) 

E 


-  !  - 

re  x  re 


*e 

■rt 

* 


c  z 

*  * 

«*4  W 

O 

Z 

IA  CO 
M 

C  4) 
4>  U 

c  re 

S.a 


X 

IA 

o; 

E 

*7 

•H 

E 


E 

a 

§ 

u 


VJ 

CO 

WJ 

CO 

VJ 

CO 

8 

X 

X 

a 

s 

u 

s 

CO 

CO 

CO 

z 

g 

g 

1 


appendix  c 


LISTING  OF  REPSL1D 


The  version  of  RPSL1D  discussed  in  this  report  is  catalogued  in  file 
RPSL1D,  cycle  2,  in  UPDATE  form  at  the  time  of  this  writing.  There  is  one 
COMDECK  called  MAIN.  This  is  listed  first.  The  remainder  of  the  listing 
is  the  COMPILE  file  image  of  RPSL1D  formed  through  UPDATE  with  the  COMDECK 
images  replaced  with  the  COMMENT  statement  *CALL  MAIN  HERE. 

This  listing  gives  the  correct  UPDATE  card  identifiers.  It  could  be 
used  for  simple  changes  through  UPDATE.  For  involved  changes,  the  user 
would  be  wise  to  use  the  FTN  compiler  to  create  symbolic  reference  maps. 
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•compfck 

COMMON  y.  n  03)  .7  (  1(13)  .TP  (103)  .0Z  (  Hi 3)  .FM  IP  (  10?)  .Fm??h  II  031  .  MAIN 

1  F**??K  (103)  .FNITM  103)  .FN?P  ( |  03)  '  FN2K  ( 1  03)  .SNH  I  In?)  .SM  (1  03)  . 

2  TEMP (1031  .P  ( 1  (13)  .FOSL  1  (103)  .EPS!  ?(l<>3)  .FPSU1  (103)  .FPSl,?(103)  .  MAIN 

3  SUM  (103. -if)  .Sit?  (1  03.3*0  .L«»T(  103.6)  .  K(U 

*  nfPSl  (N)  »  Ofc  PS?  (A)  .  ZF  T  A  ( f> )  .  ZFTASu  (f>)  .SS1MN(NI  «  SS2MN  It)  .MM  .  MAIN 

e  MCYCP(50) .NC30P (SOI .JCYNLP(Sfl)  MAIN 

COMMON  NN,  wr>22.PTD?.nFTAl  .DfcTA2.HADIUS,NN3D.DAMMF,  MAIN 

1  OP'ACT.MnAMP.TPAMP.LOAn.Pf  L6AM  MAIN 

COMMON  F  .FNO.b.PMAT.SIG2.GAM2.70.7l  .1  A  V  EH  .  PE  LT  A  T  .  T  I MF.  .LPPES*-  .  NAIF 

1  nN'N.NCYCLF.NPITF .NCONT.NSTpN.CtNfP.ClNFS.CINFP.Cl .C2.NPL0T.  MAIN 

?  PfLSG.TA.MAXf .MHITF.CA .CP.CINET.STPEN.PLAST.TNWG.  M  S.FNP.MAIN 

3  KPPINT.AnFLP.KMtSM.IMC61.IPCE2.ISP.NSFL.KJMtfX.YLnF AC.ALH  MAIN 

COMA  ON  Ml(M.Nt?(N).CM<M  «0N2  (M  .  PN  ( F  )  •  FT  AO?  (  N )  . ANGLF  (M  .  MAIN 

1  ANf-LP  (M  .NETACM6)  .FPFF1  ( M  .FPSS2  (M  .  JCPK  (3)  MAIN 

COMMON  f'N.NOl  ,N02.(,’N1  .(.N2.ni  .02.FTA02  MAIN 

COMMON  MSA (M  ,hS« tt ) ,AF6( M  .HSH (F ) .Oil  1  (M .hi?? (01 .  MAIN 

1  A1 It  <M .A??t (Ml .M  n  (Ml ,H2?1  (N) ,A1 1?(M) ,A?2? (Ml .P) 12IM .*??? (N>  MAIN 
COMMON  OSP(M) «PSm(6> .SSIGI6I  .SEPS(N) .SFIM) .SlGZSbCM  .NT ( (j  1  MAIN 

COMMON  «IG77(h).IM  MAIN 

CO“MiN  P2.72.0P2.072.P22.722.0H22.DZ2?  MAIN 

COMMON  All .A22.SPA.CS21 1 .CS222.Hll.P22.BT.PMll .SN?2 ,C1N k SI .C IN fc S?  MAIN 
COMMON  IS  MAIN 

Common  MP.N2P.N1V.N2W  MAIN 

COMMON  N1 A.N2A, IM.PT02M.H022M  MAIN 

COMMONSNpM  (103)  .SN  KM  (1(13)  ,SIMM(10?.3M  .SIP2M  ( 1  03.3M  .L»ATM(ln3.N)MAIN 
COMMON  rFFX.NnEKX.TOFFX.EPSX (H) .NFPSX ( P ) .TFPS* (HI .  MAIN 

1  »7 ( 1 03 1.ZZ( 103I.G«11 (1031 .GO?? (103) .GLl 1 ( 1 03 1 . GL 2? 1 1 03 1  MAIN 

COMMON  |iN3(M  .M3(N)  .PA1  ( 1031  .DPI  (103)  .OA2U03)  .IIP?  (103)  MAIN 

COMMON  I0AUSS.W2FTA (f ) ,»?FTSU(M  MAIN 


? 

3 

4 

5 

N 

7 

F> 

Q 

10 
n 
12 
13 
1  A 
IS 
IP 
17 
1  A 
IP 
20 
21 
22 

23 

24 

25 
2N 
27 
?H 

?o 

30 


PPuGpah  PPSL  II  (  INPUT. 01  TPl.T.PL«lT.TAPF4*  Input.TapP6*(  «»T«t, T  • 

1  TAPF3»6i_nT) 

C  *r»l_t_  PFPF 

C 

r  CPAFDrS  [N  SURFACE  STRAIN  tCNF"TAT IC6S  6/26/76 

r  OPTIONAL  ':*n.iSi>lAN  IMEohaTK-n  6/0/76.  i.Sr  ii  IK  I{-«fSS*l. 

C  t'PP“IT  IJP  TO  6  “eCHAMCAL  6(jf>L  A  YF  PS .  NSFL  .Lt.  6 

C  CI-OIPEPSIL  "ASTF.h  PBCCPAP  FOB  Cut.  10/13/77 

C  IFOFIP  FOP  PE AM  CP  SLAB  SYMMETRY.  PPtSS  Flip  CONSTANT  P»EccuBfc 

C  I“Ctl»l .?.CW3.  IflCF?al .2.3.0R4.  H/28/74  1/7/76 

T  CHANCES  TO  PFPM IT  SYPHETHY  AT  FNC  1  INSFPTED  ON  1/7/76 

f  PFPVIT  UP  TO  A  MECHANICAL  SUBLAYERS.  NSFL  .LE .  4 . 

<:  FKiP?  FPFE  IF  I  MCE2»4 .  8/2H/74.  HAY  UNL  Y  bF  liCUO  FCW  OFAMS. 

C  DIFFFPFNCES  FOP  FM  CENTEHFO  AT  HF  Sh  POINTS.  FOP  F  M  IN  MP.SPES 

C  TAPE  (NPLOT)  PLOTTING  0*TA 

C 

C  1MTIATE  PP06RAP 

C 

CO  O01  L*1.4 

fpsxu  i*-iooonn. 

ooi  EPS* (1 ♦41*100000. 

OEFXsO.O 

c 

F-PLOT*3 
PEW  I MO  3 
NLP*1 

NAN*  1 

Call  stapt 

c 

C  UELETE  PESTAPT.  II. F.  wPTAPE) 

r  IF (NCONT  .LF.  01  (iCTO  15 

C  CALL  WPTAPE  (?) 

C  NPP INT*  I  Nl  YCLF  -will!  (FCYCLEtPDELP)  I  ♦NT  ELP 

C  C-l  TO  40 

IS  F'CYCIP*0 
TIMF*0.(I 
CINFSaO.O 
rnAHP*p.o 
CINFP*0.0 
EFSaO.O 
01*0.0 
r?«o.o 


C  SFT  INITIAL  I)  I SPL ACFmEN T t PPFSbURF • STP A  IN  AN  0  STRESS  »  0 
C 

17  OC  2H  M*1,\N 
OP  (MsO.O 
I  Z  (M  *0.0 
P  CM  ■o.o 
FPSL1  INI  •«.!» 

f:psl?  ini  ■o.r* 

FPSI'l  (Man.n 
FPSn?(M»n.o 
00  2«  Ka] .KJMAJI 
SIP1  IN. A  1*0.0 
StP?<N.Kl*0.0 
SI61N (K.M  *0,0 
SI«?N{4. 01*0.0 
?«  CpNTINUF 
r 

*f  :  U-'Ai  :  -il.L 


PPCL10 

p  K  c  L  10 

TAP  , 

FrPNLlO 

PFSL  10 

PPSL  If' 

PPSLlD 

PFSL  ID 

PFSL1D 

PPSL10 

PFSLID 

ppSLin 

PPSLIO 

ppSLin 

PPSLlD 

PESLIO 

PPSLIO 

PPSLlD 

PPSLIO 

PPSLlD 

PPSLlD 

PPSLlD 

HPSL10 

PFSLID 

PPSLlD 

PPSLlD 

PPSLIO 

PFSLID 

PPSLlD 

PF  SL  in 

PPSLlD 

PFSLID 

PPSLlD 

PFSLID 

PFSLID 

PPSLlD 

PFSLID 

PPSLlD 

PPSLlD 

PPSLlD 

PPSL 10 

PPSLIO 

tsl  i  n 

PFSLID 

PFSLID 

PFSLID 

pPSI  10 

PFSLID 

PFSLID 

PF  SL  in 

PPSLlD 

PFSLID 

PFSLID 

BFSL10 

PPSLlD 

PPSLIO 

PPSLlD 

PPSLIO 

PPSLIO 

PFSLID 

vt  i  M  n 

pFSLil' 


2 

3 

1 

F 

6 

7 

P 

0 

10 

11 

1? 

13 

14 

15 

16 

17 

18 
IP 
20 
21 
2? 
?3 
?4 
25 
?6 
27 
20 
20 

30 

31 

32 

33 

34 

35 

36 

37 

38 
30 

40 

41 

42 

43 

*4 

45 

46 

47 
4H 
40 
*0 
PI 

52 

53 

54 
65 

46 

57 

SH 

6« 

60 

61 
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n  n  n 


ntLG*“«rf 
CALL  iKiFoh 
r all  sthain 

Lj{.  31  .'■■M  r « .\  ?P 

a??  *  p.5*ir»?iM*r»?(A-n) 

IE(h  .FO.  h» Nt A??*C A  5 <M 

IF(h  .FfJ.  N2H)  A2?»0A?<N-n 

FHU  (M  *  1 .0/ lf.«l  (M-?.0*7U*DP1  CM  1 

<H'?2(E>  *  1.0/ (  A22-?.04ZO«Ort2 <N)  1 

gluim  *  i .o/ir.Ai  im-2.o*zi  •ruiiM  t 

C-1.221N)  *  1.0/1  A?2-?.P*ZI  *0H?(M  ) 

I.A11M)  *  0.0 

oa?in-ii  »  o.n 

DPI (N)  *  0.0 
DH2<M  s  O.o 
33  COMIM'F 

CA?<N2P>  «  0.0 
f.  FSU?  CN?Ei )  =  O.o 
E PSL2 <N?F  1  *  0.0 
IF1LOAP  .OT.  OlfiOTl)  43 
f 

CALL  InVFL 
C 

nil  34  N»MV.N2V 
cp  (mi  »ofltat*i;p  ini 
P7  (N ) «0ELT  AT*P7  <h) 

35  CONTINUE 

PALL  POONCI 
CALL  K1AFT 
CINFS*2.0»CINFT 
TNRC-»CINPS 
C 

IF (LOAD!  42, *5. 43 
C 

4?  CALL  PP'OPF 
C 

43  PO  44  MMV<4?« 

I.H  (N)  «f>R  <N)-P  (E  )  *SNM  (Ml *T€NP (N 1 
CZ  <N) »OZ (N) -P  (M  *SNK  (M  ‘TENPIN) 

44  CONTINUE 
CALL  POONMJ 
CALL  PWOma 
«N  9*FM« 

CALL  FINET 
TNRG»CINET 

4*  IFINCYCHlll  .Fly.  0)NNN«? 

C  NPITF  INITIAL  CARTESIAN  COORDINATES.  PRFSSUHF 

ROITF  ( (•  •  3 (> 0 ) 

4HTF1A.4I  0)  (N.P(N).7(M.P(M«N«1.NM 
C 

4«  call  foata  m 

END  initialization 

so  NCVClE«NCYCLt»l 
T I NF  »T 1 NE *OfLTAT 
C  CPFCF  EON  final  step 

IFCWVCIF  .OT.  NA*C)C,(.TO  70 
C  ChFCF  IF  CAtL  PRESS  IS  NFFC'Ffl 

IF1LPOFSS  .PE.  NCYCLE1CALL  PPFSS 
CALL  POSITN 
■ 1  l  '  ’  •• 

CAIL  STRAIN 


*P?L10  (-4 

RPSLID  n? 

PPSLIO  An 

-i-SLir.  *-7 

HPSLIO  88 

PE  SLID  89 

PPSLIO  70 

RESL1U  71 

PPSL10  72 

NPSL10  73 

REEL  10  74 

RFSLIP  75 

PPSLIO  7 t 

PPSLIO  77 

Pt-SLIU  7E* 

HP  SL 10  TR 

HPSL10  RO 

RE  SL 1 0  HI 

HPSl  10  82 

RPSLID  83 

HP  SL  1 0  84 

PPSL10  85 

PE  ELIO  R8 

HESL10  87 

NPSL1P  88 

PPSLIO  MR 

PFSL1U  SO 

NPSLIO  SI 

RPSLID  42 

hPSLIO  93 

PPSLIO  94 

HPSLIO  95 

RPSLID  98 

PPSLIO  97 

HPSLIO  98 

HPSL1D  99 

HfSLIO  100 

RPSLID  101 

hPSLIO  102 

HFSLIO  103 

PPSLIO  104 

RPSLID  105 

RPSLID  IH8 

RPSLll)  107 

RPSL10  108 

RPSLID  109 

HPSLIO  HO 

hPSLIO  111 

hPSLIO  1)2 

HPSLIO  113 

HPSLIO  114 

PPSLIO  HR 

NPSLIO  lit 

HPSLIO  117 

PPSLIO  HP 

PPSLIO  114 

PPSLIO  120 

HPSLIO  121 

PPSLIO  12? 

PPSLIO  1?3 

PPSLIO 

-  r  w  !  .  *'•  ;  >  ~ 

hPSLIO  1/8 
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CALL  MOTION 
CALL  P  0  A  T  A  (?) 

CALI,  DAMP 

r  CF>Ft>  POP  RESTART  (  i-lap 

C  IF  (NCYCLF  .ml.  nPITF.  )I-«'.T0  Si) 

r  CALL  APT  ARE  (1) 

C  «FITF*.vPITt*KfITK 

r  call  pcata  (3) 

goto  so 

70  IE (NCYCLfc  .LT.  ?)5TCP«  MAIN.  NCYCLE  .LI, 

CALL  PC ATA  O) 

CALL  POATA  (4) 

STOP*  PPOGPAR  COpPLETF* 

300  FORMAT  ( |  hi  .?0X  .1«HIMT!AL  COORDINATf  S* ?3 X  * PHPRFSSUPF // 
l  3X. ]HN. 1 1X.4RR (N) ,?] * .AN?  (N ) ,?1X  ,4PP  (N ) ) 

400  FOPPAT ( ( 14.3 (?X.K?3.1M ) ) 

FNO 

SI'PWOUTINF  START 

C  *C  ALL  PAIN  HERE 

niPFNSION  TITLE (H) 

PE  AO ( c • | Oft )  TITLE 
*>R  ITF  (6. 140)  TITLE 

FEAn (S. | OS) NOOSE .NPFSH.lAYEP.YLOFAC 
WFAO(S. 10S)  pAXC.nCOnT.FRITE.OELTAT 

PEAms.uo)  noose  . ipce?.nohsf.ipcei 

PF  AC)  ( S  .  |  1SI  LOAO.LPPESS.PCAPF  .DAPPF  .['FACT 
PF  AO (S , 1 ?0 )  E .FMl.S IfiZ  «PPO*TPICKN.n<FL  « ISM 
NP’ITF  ( A  *  1  NO )  E  .FnU,SIG7«PH0.TEiICKN 
RP  I TF <  F  «  1 70 )  NCONT.RAXC .NPWINT.NPITE 
PF  I TF  (  A.  1  7s)  L  A  YEP  .NSTPN.LOAII.LPPFSS 
WPITF(A.IPO)  IRCE  1 *  IRCE? 

IE  (NSFL  .tO.  1  •  ANT .  ISP  .FU.  0 )  GOTO  700 
IE (NSFL  .EO.  1)  IfiOTO  700 

PtAP  ( 5. 1 ?S)  (SSUMb) .SFRS(d) *OSP(d> *HSP(d> *d«l.NSFL> 

700  IF (NSFL  .LT.  1)  ISRa-l 
IF(NSE'l  .IT.  1)  NSFL*1 
PF  »0(S.  HO)  NPPINT.  ( JCHX  ( d)  . d»l  .3) 

PE AO(S.) 10)  NLPCY, (NCYCP (d) *d«l.MtPCY> 

READ  15. 1 10)  NLPWIN.  (dCYNLP(d)  «d«l  .M.PPIN) 

PFAfMS.1  10)  N3P.  (NC30P  (  d)  .d»l.N30> 

SSItt(l)»sir,7 
SF  DS ( 1 ) «S 1 07/F 
A  d  P  A  X  * I  oYEP«n«FL 
sr  i  :  i  =E 

SIG7SCJ  ( 1  )  «SSTG  ( 1  )•*? 

SIPZ7C1)  ■  SSK-U) 

TO  70S  jbI.nSFL 
if (ISP  .LT.  I  HiOTP  7CA 

IF  ( PSF  ( d )  .07  ,  o.o  ,AN(  .  PSP(d)  .NT.  0.0)  GOT  l'  7»3 

7R1  rrITF (N.7o?) 

70?  F  PRMAT  ( //*A)H  F  F  POP  IN  STRAIN  HAPPFFlM*  ON  STRAIN  NATE  );A  T  A 
STOP*  START.  ERROR  IN  STRAIN  DATA* 

703  RS«(d>  «  I.O/PSR(d) 

704  IF (d  .FO.  1)  Gt To  70S 

IF (SFRS(d) .LF ,SEPS(o-l ) )  GOTO  791 
SF (d)*(SSfft(»)-SStG (d”l ) ) / (SEES (d)-SERS (d-l > ) 

NT (d-I ) «<SE (d-l ) -SF 
SI07SL (d)«(€*SFPS(o) )*»?  ' 

SIG77(d)*(F*SEPS(d) ) 

70s  CONTINUE 

RTINSFl )*SF(FSFL)/F 
r  set  i‘p  bni  NOnUY  nintooi  s 

i*  iRCEi®l.?."Pi.  .Fl,rP*i»c...  *  t.  -  **/  r  R  /  V  ■*  1  /  l  /  i  h 


rREI.  Ill 

l  ?  7 

PPSL10 

I  ?R 

PF  EL  1  n 

1  ?u 

pf  el  1 0 

130 

PFSL  10 

131 

rfelio 

13? 

hF  EL  1 1 1 

133 

RFSL10 

134 

pfslid 

135 

pfslio 

13F 

WFSLIO 

137 

rfslid 

13H 

PFSLID 

130 

PFSLIO 

140 

PPSLID 

1*1 

PFSLID 

14? 

PESLID 

143 

START 

? 

TAF 

1 

START 

4 

START 

4 

STArT 

N 

STA  kT 

7 

start 

P 

STaRT 

9 

START 

10 

ST/PT 

11 

START 

1? 

START 

13 

START 

14 

START 

15 

START 

IN 

START 

17 

START 

1 P 

START 

19 

START 

20 

start 

21 

start 

2? 

START 

23 

START 

2* 

START  ' 

25 

START 

?N 

start 

27 

»T  FR  1 

r  b 

staft 

?o 

START 

30 

START 

31 

START 

3? 

START 

33 

start 

34 

STArT 

35 

start 

3N 

start 

37 

start 

3R 

START 

30 

START 

40 

start 

*1 

STArT 

4? 

STArT 

43 

STArT 

44 

start 

44 

<-Tf  FT 

.1  i  4  r  r 

-  / 

75 


n  n 


B 


h 


* 


4 


l 


FT  AFT 

4H 

f?PafiP7FR«..ilr 

START 

49 

A  1  VaR)  1 

PTART 

*  —  ( 

ciu-T 

4) 

K  1  A  a  (v  )  v  - 1 

ST/rT 

F? 

7  ?  Aa*.  *  1 

ST  APT 

S3 

I  f  (  (POf  1  .  ►)  .  ,>)  |\  1  vs  v  |  R 

START 

Pp 

17  ( IFI  f  ?  ,F  <:.?  .OP  ,  I  ACE  ?.►  (■  .4  )  MPI/aK^R 

FT/  h  T 

FF 

7-f  ah?A*) 

STAPT 

FA 

IF ( IR( F?  .FF.  4)  )>KaM?P 

STArT 

57 

PFAMF.IJo)  SCOSfc .FTAO?.NSTRN 

start 

5P 

RF»(>  (5.  15)  (NdJSF .FTA62  (  I  )  . ANfil F  (  I  ) 

•A66LP  (  I  )  .NFTAG  ( I  )  •  I a  1  •  NSTP6' ) 

STAPT 

59 

CALL  I  NO-FOP 

STAPT 

60 

CaII  INFORM 

STApT 

M 

SIPFAIS  TC  RP06rAp  FRO*  IN«EOr. 

START 

62 

AMAl  SYRRFTPV  (RADIUS  >  0.0, 

iRall). 

START 

6.3 

FI  Ap  SYMMETRY  ( P  AGIOS  a  O.li 

IP*0  )  •  PEAR  (lfi>0). 

STAPT 

64 

IS  a  1 

START 

65 

IF  (IF  .C-T.  0)  RADIOS  a  0.0 

STAPT 

66 

IF (RADIUS  .LF.  0.0)  T5a? 

STAPT 

67 

IF  (  IF)  .fiT.  OlDFTAlal.O 

STAPT 

66 

t-NaFLOAT  (MH)  ♦HAOF/PFTA? 

. 

START 

69 

»  C  I  a(.lf 

STAPT 

70 

TP?  at  (J1  ♦  1 

START 

71 

CMaUM-KLOAT  (M,  1  1 

START 

7? 

OP  ?aFI  OAT  (Mi)?  )  -I  N 

STAPT 

73 

fC  ?0  lal.NSTBN 

START 

74 

R7  (I)aFI  OAT(MR)  *E  t  AG?  (  I )  /OK  T  A? 

STAPT 

75 

A 11  (I)apf  (I) 

START 

76 

A' !?(!)■►  Il  <  t )  *  1 

START 

77 

FH«M1IH 

STAPT 

7« 

P7?aMT  ? ( I ) 

START 

79 

I’M  (  I )  aMi  (  I )  .PM) 

STAPT 

»0 

P7?( I)aPF?-PK ( f) 

START 

PI 

IF  (ON  1(1)  ,GT .  O.SK-CTO  16 

START 

P? 

OF  1  { I )  LF  1/?. 

STAPT 

63 

Midi  s  Ml  <T)-1 

START 

64 

063(1)  a  f»M  (I)  ♦  n.s 

STAPT 

«5 

l-OTO  17 

STaPT 

H6 

OA'l  (I)  <=T  )/? 

start 

67 

)6 

7  13(1)  a  Hill) 

start 

8P 

IF  3(1)  «  Mllll  -  O.F 

START 

»9 

1  7 

IF  (Midi  .r,T.  1)1, r, TO  in 

START 

90 

A  I  3  (  n  L  *  ?.  >R  ,  ft  TO  ?. 

S  1  ^  p  i 

-1 

Midi  »  ? 

START 

9? 

rum  a  o ,n 

start 

P3 

or  to  ?o 

STAPT 

94 

1  R 

IF(MJIT)  .  L  T  .  r-rTo  20 

STAPT 

94 

M  3  (  I  )  r,T  6?R-?.  Shift  TO  N?R-?. 

START 

06 

7  13(1)  a  ('<?«-? 

start 

97 

rifim  »  i.o 

STAPT 

9F 

?0 

LONTlKhF 

START 

Op 

7AYFPa(  AV7R 

STAPT 

10  0 

“UCOSAM  CONSTANTS 

start 

10) 

IF  (  l  H  .16.  0)  /l.'a  l,  .s  »Th  1  CAN 

START 

10? 

7t  *7U»Th  I  (  KM 

STAPT 

103 

<-AM/aTH  Iff  M»M-n 

START 

104 

7  TFoPaNPR  TNT 

start 

1  OF 

PP  ITF  aMf(.M  T *7 R  I  TF 

start 

1  06 

S0“*P«O.0 

STApT 

107 

SI'P-?ABa(i.o 

START 

10H 

O-JiAFit  TC  C|ltM||  tr.  v  i;,M,CC(  Ik  >«■■' 

TFi-UATroc  i>pl  7rTONFFS  . 

ct  /  PT 

’  OP 

i7  i  /  s  ,r  c.  ,  |  i  1  •!» 

>  I  f  l>  ' 

1  1  0 

4 


76 


i 


IF  <  1  »'Al  ««*  ,NF.  norTO  R(  START 

St  T  u(-  FOR  GAUSS  [  An  I  R  T  t  <?»•  AT  t  UR  »Ilh  lAYHh  (.(_£.*■>  K  1NTS  ST/hT 

PCTfUAl  .4?.4?.44.AC,4A|  .|_AYFR  START 

Al  7fc  T  A >11  X  n.r  RTAhT 

Mil  :  M'  S  T  A  *  T 

GOTO  SO  START 

4?  ZFTAl?)  x  1.-77‘’S.i;tSl-Lf.£h«7l,  ST.RT 

7F  T 4  (  t  )  X  -7M.IZI  START 

k(l>  x  1.0  ST  A  ST 

►  (?>  x  1.1)  START 

GOTO  SO  STaRT 

a?  7FTA(3>  «  0.774‘>‘»<‘>ARO24l4fi3*7ll  START 

Z£TA (2)  s  0.0  STahT 

7F  T  A  { 1  )  x  -2£TA(3>  START 

V(l)  x  O.SRSRRRSSSRCRRSH  START 

«(?)  x  n.SMRFHMHHrtHHAHHO  START 

►(3)  x  S(l!  START 

GOTO  SO  START 

A*  Z£Ta ( A  1  x  0  .HA  1  1 3i>?  1 1  HSaOS3*7d  START 

7FT4I3)  x  0.33SSi«lflA3SPA856*7U  START 

7FTAI?)  x  -ZFTA(3)  START 

7F Ta ( 1 ) x  -ZKTA(A)  -  START 

*(t)  *  0.347HSaRaS1 37484  START 

W(?)  s  o.6SZ1ASi^*ha?s<.8  START 

w(3>  x  *  (?)  START 

M4I  «  Mil  START 

or  T<i  SO  START 

AS  7FTA  (S)  c  0.RDM7SR4SR3HR6AA7II  START 

ZFTA(4)  x  0.S3HARR31 010SGH3*ZU  START 

7FTA(31  x  0.0  START 

ZFTAl?)  x  -ZfcTA(A)  START 

ZFTA(l)  x  -ZETA(S)  START 

R<1)  X  A.?1RRi>HHHR0SMRS  START 

M?)  x  0.47SR?8h70Aqo3Hf  START 

» ( 3 )  x  O.SSHRHOhRRHHHHHR  START 

A ( A )  x  k(?)  START 

R (S)  x  k ( 1)  STAhT 

GOTO  RO  START 

aa  2FTA(M  X  (I.S32arRS1a?031S?*Zii  START 

7F  T  A  ( S  )  x  O.GRl?()R?H6AAR?5R*ZO  START 

7FT«(a)  x  o.?3RAI«1Raor3 J«7*7i,  START 

7t  T A ( 3 )  x  -7F  T  A ( A )  START 

7FTA (?)  x  -ZFTA(S)  START 
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APPENDIX  D 


PLOTTING  PROGRAM  FOR  RPSL1D 


The  plotting  program  is  an  independent  program.  The  RPSL1D  program 
stores  data  for  plotting  through  unit  3  which  is  equivalenced  to  file  BLOT 
on  the  PROGRAM  card.  The  plotting  program,  which  we  have  catalogued  in  file 
RPSL1DPL0T,  reads  data  from  unit  3  also  equivalenced  to  file  BLOT.  For  some 
runs  it  may  be  advantageous  to  catalogue  BLOT  and  repeat  the  plotting  program 
to  reproduce  plots  that  were  unsatisfactory  because  of  scaling  (input  card) 
or  system  problems. 

This  program  uses  the  BRL  plotting  subroutines  described  in  ARDC  TR6 
(BRLESC  FORTRAN  PLOTTING  SUBROUTINES,  Monte  W.  Coleman,  John  V.  Lanahan, 

July  1970)  as  amended  for  CDC  by  local  publication  SPB-6-78,  May  2,  1978. 
Conversion  to  SCOOP,  the  plotting  system  used  in  Reference  2,  or  another 
plotting  system,  would  not  be  difficult. 

The  main  program,  RP1PLT,  reads  binary  data  from  unit  3  (equated  to 
BLOT)  and  controls  the  program  flow.  If  the  input  variable  IFLAG  is  1,  an 
array  of  data  for  one  time  point  is  read  and  stored.  If  IFLAG  equals  2, 
different  data  is  read  and  PL0T3D  is  called.  If  IFLAG  equals  99999,  or 
certain  abnormal  conditions  occur,  the  program  calls  subroutine  GRAPH. 

Subroutine  PLOT3D  reads  a  control  card  on  the  initial  entry  containing 
DEFLM,  SOFC,  and  SF  with  FORMAT  (3E10.3),  and  sets  up  scaling  for  isometric 
plots  of  the  "center  line  deflection  profile"  of  Reference  2.  This  is  a 
plot  of  the  initial  position  of  the  reference  curve,  and  the  current  position 
of  the  curve  with  the  difference  magnified  by  DEFLM.  If  SF  is  greater  than 
zero,  it  is  the  reciprocal  of  the  scale  factor  between  the  internal  length 
measure  and  inches  on  the  plotting  surface.  If  the  input  SF  is  not  greater 
than  zero,  the  program  assigns  an  SF  that  attempts  to  scale  the  plot  of  the 
initial  curve  into  an  SOFC  by  SOFC  inch  square  which  also  includes  the  origin. 

Subroutine  GRAPH  plots  displacement  increments  DR  and  DZ  vs.  time  at 
point  (ETAD1,  ETAD2)  and  the  energy  balance  plot:  time  vs.  kinetic  energy, 
kinetic  energy  plus  strain  energy,  total  energy,  and  total  energy  less  damping 
work.  Subroutine  GRAPH  calls  subroutine  STRAIN  to  produce  NSTRN  plots  of 
strains  in  the  coordinate  directions  vs.  time  at  points  prescribed  on  input 
cards  13.  If  the  PLOTP  option  is  included,  GRAPH  also  produces  the  NNPE 
prescribed  plots  of  P(N)  vs.  time. 

All  the  COMMON  variables  are  included  in  COMDECK  MAIN.  The  longer  arrays 
are  put  in  LEVEL  2  (a  special  version  of  FIXSCA  was  inserted  to  use  them);  if 
more  than  3000  time  cycles  are  recorded,  these  arrays  and  MAXC  must  be 
increased  to  plot  them. 

The  following  list  is  the  COMPILE  file  image  of  RPSL1DPL0T  formed  through 
UPDATE.  This  listing  gives  the  correct  UPDATE  card  identifiers. 


r>  r»  n  n  r» 


FPOGPAP  PHPLT (INPLT.OLTPUT.HLOT.TAPt 1 3  .  T AHF5* l NPOT . T APF 8*0UTPUT . 

pPJPLT 

2 

1  TAPF?«ei'Vn 

PFJpLT 

3 

PLOTTING  FOP  1-0  REFS IL . 

PH  PI  T 

a 

PFPSIO  PLOT T INK  PACKAGE  (CALCOPP  PLOTTING. 

SPF-G-7P) 

PFJPLT 

5 

IF  Nl:YCL>PA*C.  INCREASE  MAXC  AM)  APPAYS  IN  PARPAY  AN  D  PLOTP. 

PPJPLT 

8 

PLOTTING  PIN).  9/30/75  (pn(3ou?.9)> 

PP 1  PL  T 

7 

pfiplt 

8 

COMMON  E  TAG? (8 ) ,PN<8)  .NETAG IP)  .NSTPN.HAXC 

PAIN 

2 

COUPON  Bp ( 1 03 1 .72(103) .Y 1 103) .7(103) 

PAIN 

3 

COMMON  OAT (20) .NCYCLE.TIPE.  ETA02. 

CIN. NCYCL 

PAIN 

A 

COMMON  NM>E.NPE(91  ,Pf)AT<91 

PAIN 

8 

COMMON  /PAPSAY/ 

PAIN 

8 

1  TIPI  300?).  U2(  3002).  03 1  30  02).  CINI 

3002),  STC (  3002), 

PAIN 

7 

2  TNP<  3002) .OAPPLT (  3002 1 .SPSS 1 U PO 12) .EPSS2 (180 12) 

PAIN 

8 

CCMPON/PLOTP/PD (  300?, 0) 

PAIN 

9 

LEVEL2.TIM.PC 

PAIN 

10 

MAXC  «  300? 

PPIPLT 

10 

NNPF  a  0 

PP1PLT 

n 

NCYCL*0 

PPIPLT 

12 

NPL0T*3 

PP1PI.T 

13 

BEHIND  NPLOT 

RP1PLT 

14 

PPIPLT 

15 

PFAO(NPIOT)  ETA02.0N.NSTRN 

PHPLT 

18 

BEAD (NPLOT)  (ETAG2II) ,PM I ) .NET AG  1 1 ) ,1*1 .NSTHN) 

PPIPLT 

17 

RF AO (NPLOT INNPE. (NPFII) , I *1 .NNPE )  ACTIVATE 

FQH  PLOTP 

PPIPLT 

18 

Of  All  (NPLOT)  NCYCLE.TIPF . N 1 . < PR (N ) , 7? IN) .Na 1 . N 1) 

PPIPLT 

19 

II*2#NSTPP*« 

PPIPLT 

20 

PPJPLT 

21 

SAVE  INITIAL  SHAPE  NEEDED  FOP  OFFLECTICN  MAGNIFICATION  IN  PI  f.T3D 

PPJPLT 

22 

DO  5  Nal.M 

PPIPLT 

23 

YIN) *72  IN) 

PHPLT 

24 

7  IN) *PP (N) 

PPJPLT 

25 

5  CONTINUE 

PHPLT 

26 

CALL  PL0T30IM) 

PHPLT 

27 

RP1PLT 

28 

10  BEAOINPLUT)  iflag 

PPJPLT 

29 

IF (EOF (NPLOT)  .NE.  0)  GOTO  28 

PPIPLT 

30 

IF  1 IFLAG  ,E0.  999991  GOTO  30 

PPIPLT 

31 

IFIIFLAP  ,EQ,  1 ) GOT  20 

PHPLT 

3? 

IF  1 IFLAG  .EO.  2) GOTO  25 

PPJPLT 

33 

HPITFI6.9H)  IFLAG, NCYCL.NCYCLF 

PPJPLT 

34 

911  FOPNATI//'  PAP  SIGNAL  FROM  TAPE.  IFLAG .NCYCL ,NC YCLE * • , 3 11 0 ) 

PPIPLT 

35 

IF  PUN  rAlLa,  GF  T  PARTIAL  PLOT. 

PP i PL  T 

34 

GOTO  28 

PPJPLT 

37 

PPJPLT 

38 

?fl  ME  AO  (NPIOT)  NCYCLF ,(DAT(1),I*1,I1) 

PPJPLT 

39 

1  . (POAT(d) .J*l .NNPF )  ACTIVATF 

FUH  PLOTP 

PPJPLT 

»0 

IFIFOF(NKOT)  ,NE.  0)  GOTO  28 

PPIPLT 

41 

NCVCL»NCYCL ♦ 1 

PHPLT 

4? 

IF  (NCYCL  .GT.  NAAC1G0TC  2» 

PHPLT 

43 

TIM (NCYCL 1 *0A  T 1 1 1 

PPIPLT 

4* 

1)3  I  NCYCL  )  *C  AT  (?) 

PHPLT 

45 

I'2(NCYCL)  «UAT  (3) 

RF  JPLT 

48 

CIN  (NCYCL ) *PAT 15) 

PPJPLT 

47 

STC  (NCYCL  )  *f!AT  (6) 

PPIPLT 

48 

TNP  (NCYCL ) *CAT (7) 

PPJPLT 

49 

OANPl.T  (NCYCL )  "OAT  ( 8  ) 

PHPLT 

50 

DO  2?  1*9.11.2 

PPJPLT 

51 

vJ*NCYCL  *  n»XC«<  1-9) /2 

RPJPL  T 

52 

FPSSI (d)«OAT(  II 

RHPLT 

53 

fuss?'.m*D»7?!*'  ) 

UP  1  PI  T 

84 

22  CONTINUE 

pmpl  r 

-8 

98 


IFINNPF  ,fo.  OIGOTO  ip 

PF1FLT 

ro  23  I»1  .nnPE 

RF  1PLT 

N7 

?3  ppincycl.ii  »  pdatui 

WHIPLT 

SM 

Pete  io 

PF1FLT 

2*  PE  AO  (NPLOT  )  NCYCLF .T IMf ,N1 . (Wp (N) .77  IN) .N*l .Nl ) 

PF1PLT 

NO 

IF (EOF (NPLOT)  .ME.  P)  iPTO  2F 

RFIFI T 

Ml 

CALL  PLOT3D  (M  ) 

PUPLT 

GOTO  IP 

PF 1PLT 

N3 

ppiplt 

64 

?h  HR ITE  <  6 , 1 00 )  MAXC 

ppiplt 

NS 

IPO  FORMAT I//20H  FPROH  NCYCL  >  MAXC«.I5/) 

RP1PLT 

6H 

THE  LAST  SET  OF  0»TA  IS  NO  GOOD.  PICT  THE  WEST. 

PPIPLT 

67 

NCYCL  ■  NCYCL  -  1 

PPIPLT 

6H 

ppiplt 

69 

30  CALL  ghaph 

PPIPLT 

70 

CALL  PLTPGE 

PUPLT 

71 

TAX  EXIT 

RF 1PLT 

72 

END 

PF1PLT 

73 

SUBROUTINE  PL0T3U (I?) 

PL0T3D 

2 

COMMON  f  TAG2 (6 ) .PN<6) «Nf  T AG  (6 1 .NSTPN.MAXC 

MAIN 

2 

COMMON  RP< 1031 ,77(103) *Y (103) .Z (103) 

MAIN 

3 

COMMON  OAT ( 20 ) .NCYCLE « TIME.  ETAD2.  ON, NCYCL 

MAIN 

4 

COMMON  nNPE.NPE(W)  .PflATIS) 

MAIN 

5 

COMMON  /PARPAY/ 

MAIN 

6 

1  TIM<  3002).  U2<  3002) .  U 3  (  3002).  CIN(  3002).  STC< 

3002).  MAIN 

7 

2  TNR<  3002) .OAMPLT <  3002) .EPSS1 (10012) .EPSS2II00I2) 

MAIN 

H 

COMMON/PI.OTP/PO  <  30  02.9) 

MAIN 

Q 

LEVEL  2 • T IM .PO 

MAIN 

DIMENSION  X) <103).X?<103).  HE A 02 (3) *HEAD3(2) 

PI 0T3D 

4 

P I MENS  I ON  LAME L  <  4 ) 

PI  CT30 

C 

DATA  LAPF.L/*  wCPTMAN  •  , *P309 • « • X3979  * • »HP 1PLT  •  / 

PI 0T30 

6 

OATA  1/0/ 

FL0T30 

7 

DATA  < HF A02(IK) , IKsl.31/1 OHOEF LECTION. 10H  MAGMF IEC . 1 H>/ 

PLOT 3D 

A 

DATA  (HFAP3 (IK) . IKs 1 , 2 1 / 1  OHM ICROSt CON. 3H0S>/ 

PLOT3D 

9 

IF<I.F«.U)GOTO  10 

PL  CT3P 

10 

I«I*1 

PL<  T3D 

11 

YRAPNsYPAFNUO.O 

PLOT3D 

12 

IE  < I .LF .3) GOTC  30 

PL0T3D 

13 

CALL  PLTPPF 

PL  <  T3D 

14 

GOTO  20 

PL  T  3D 

IS 

10  XP AOF *12.0 

PLCT3D 

16 

HEAD (S,  1  1 )  ntFLH.SOFC.SF 

PL0T30 

17 

11  FOPMaT<3E10.3,Ik) 

PI  (  T30 

VMAXaYU) 

Pll T30 

20 

YMIN.YU  ) 

PLCT3D 

21 

ZwAX«7  < 1 1 

PLCT30 

22 

?M IN*7  <  1  I 

PL  (  T3D 

23 

f.O  12  F  *1  .  12 

FI 0T30 

24 

VMAX«AMAAl<V(M.VMAX) 

PL<  T30 

?S 

YH  I  Ns  A  M  I  N  1  <Y(M  .  YM  I  N  > 

PL0T3I) 

26 

7MAXSAMAX1  <7<N) , 7MA  X ) 

PI  (  T3D 

27 

7H INsAMINl  (7  <N ) »  2M IN ) 

P10T3D 

2  P 

1?  CONTINtlF 

PH  T3D 

29 

IF (SE  .NF.  0.0) GOTO  13 

PI  C  T3D 

3f 

YSs(YMAX-YMIN) /SOFT 

PI OT3D 

3) 

7S»(7"AX-7MIN) /SOFT 

PUT  30 

32 

SFsAMAXl (  YS.7S) 

PLOT 30 

3? 

IFISF  ,GT.  0.0  .*i*0.  SF  ,LT.  0.<UJ)GOT()  13 

PL l T30 

34 

SFsA INT ( SF ) 

PLCT3D 

3s 

IF<SF  .LT.  1.0)  SFsl.0 

PI CT30 

3* 

'  “I  *M.j.  3'  -u,  .?e,.’  .  ’.n.i  .‘.pt:  ' 

“l  ;  ’  ■r' 

:  7 

CALL  Pt  TSCA  <  1 .0. 1 .0 .0.0. 0.0. 1 .0.1 .0) 

PI CT30 

3P 

>9 


CAM.  PLTSYM  (n.  1  ,HF  AO?  ( 1  )  .  0.0.  0.(1. -.3) 

EFCOPF  <12.  15. STEP)  PFFLP 
15  F0«WAT(F 10.3.2H  ») 

CALL  PLTSY«  <0.1. STEP. 0.1.  2.1.-.J) 

CALL  PLTSYM  <0.1  .PFSCAl  F  l/>.0.0.  0.0. -.5) 

IFFsIF'TISF  ) 

EFCOI  F  ( A  .  1 h.STFF)  IFF 
1*  FCPPAT < I?.?H  >) 

CALL  PLTSY*  < 0 . 1  . STEP . 0 . 0 ,  0.0. -.51 

SF*1 .O/SF 

?0  YFAPN«0.0 
1*1 

Yfi  YHARs*  .O.YRARF 
XMR*1 .0 

-  CFKTEH-LIMF  DFFLECT 10F  PROFILE  - — — 

YMC*YHAP 

CALL  PLTSCA (XPAR.YHAC.O.O.O.O. 1.0. 1 .0) 

XI  I  l  I *0.0 
X?  <  1 ) *-SF*7  < 1 ) 

XI  (?) =0.0 
X? <?)*SF*Z < 1 ) 

CALL  PLTOTS < 1.0.X1 < 1 ) ,X?< 1 ) .2.0) 

x2<n*n.o 

XI  <2)*SF*Y (12) 

X2 ( ? ) *0 . 0 

CALL  PLTPTSd.O.Xl  U).X?<1).?.0) 

k*o 

DO  AO  N*l.t? 

K*K*  1 

XI  <K)*SF*Y(N) 

X?  <K) *SF»Z IN) 

AO  CONTINUE 

CALL  PLTPTSIA.O.XI <l).X?(l).K.O) 

CALL  PLTDTSO.O.Xl  Cl)  ,X?(1)  .0.0) 
k*o 

CO  70  N*1.I? 

F*K»1 

XI  (K)*SF*<Y  <N)  ♦OEFl**<77  <N)-Y  <N)  )) 

X2(K1«SF*(7(N)  ♦DEFIP*(PMN)-7(N>  )  ) 

70  fOMIhOE 

CALL  PLT0TS<2.0,X1 <1 ) ,»?<1) .5.0) 

•  •••••  •  •  •  *  •  •  «  •  »  •  • 
CALL  PLTSY«*  l.l.AHC  YClf>,0.0.  .2. -3.0; 

FNCOPF  <0.101, STEP)  NCYCLE 
101  FCRPAT ( IS. 1H>) 

XPITF ( A , 1  OR ) NCYCLE 

] OP  FCPPAT ( •  ISOPFThIC  PLOT.  CYCLE  *  1 A ) 

CALL  PLTSYM  <  .  1  . STEF  . 0 . 0 .  .A, -3.0) 

TIaTIMFM  .OFf 
EF  COPE  <  10.105.STFP)  Tl 
105  FCMHATCFH.1.8F-  >) 

CALL  PLTSYP  (  .  1  .STEF, 0.0,  2) 

CALL  PLTSYP  (.1, HEAP 3(1) ,0.0.  1.1, -3.2) 

******  •  •  •  •  •  *  •  •  *  •  • 

RETURN 

FNO 

SLRROOTIF.F  GRAPH 

CONPON  f  TAG?  <  A  )  .PN(A)  ,  At  TAG  <  A  )  ,NSTPF«PAXC 
CCMPON  PP<10?) .72(103) , Y  < 1 0 3 ) . Z ( 1 0 3 ) 

COPPPN  OAT  <?0)  .NCYCLF  .TIPt  .  ETAD?.  GN.FCYtl 

rpy.iiF  |'j)  .F^jT  'qi 

COPPON  /PaROAY/ 


FI 0T3D 
PLLT3D 
FI  (  T 30 
c)  0T3C 
Pt  C  T  3P 
H  CT3I) 
PL  < T3P 
Rt  0T30 
PICT  31) 
Pl( T30 
PL0T3D 
PI PT30 
Pl( T3D 
Pt  0T30 
PLCT30 
PLCT3P 
PL0T3D 
PI 0T3D 
PIPT3D 
PICT 3D 
PLCT3D 
PI OT  30 
PI  <  T3D 
PKT3D 
PL0T3D 
Pit  T3D 
PLCT30 
PLCT3D 
PUT30 
PI ( T3D 
PI 0T3D 
PI  0T30 
PLl T3D 
PI CT30 
PLPT30 
PI PT3D 
PL<  T30 
PLCT3D 
PICT30 
PLCT3D 
PI  f  T3D 
PI  C  T3D 
Pt  CT3D 
PL0T3L 
PH  T30 
PI 0T30 
PI CT30 
P|  f  T  30 
PI  (  T3C 
PI  t  T3P 
PI PT3D 
PI  0T3D 
Pit  T 31' 
PLl T3D 
PI t  T3D 
PI  t  T3P 
PLPT3P 
GRAPH 
PAIN 
NAIF 
PA  1  F 
v  ft  T  N- 
PA  IF 


1  tim  3po?)  •  u?(  30021.  l3  (  30021.  cim  300?).  ST<:<  3g*>2). 

2  TEB  (  3002) .GAPPLT <  30 02 ) . fcPSS 1 (1  ft 0 1 2 ) .EPSS2 (1  HO  1 2 ) 

CCMMGN/Pi  CTP/PO <  300P. C) 

LFVFL2.TIH.oc 

CImFHSION  SVMl(A)  .SYM2(?>  .  FTA?<3> .SY*<3> 

DATA  (SYHl  ( I ) . 1  =  1 ,4) /10HC0MP0ME NT  • 1 OHOF  VkCTOW  .  I OHT  I  SHI ACFME  »  3hNT 
!>/.  CSVh?(I)  .  1=1.2) /lOHTIMt  <EICB.1CH0SEC0M>S)>/ 

PATA (FTA?( n . 1=1 *3) /10HFTA?  =  ,10b  (  .SH  )>/• 

2  ( SY* ( I )  «I  =  1.2)/10H  FNFEGY.OH  PALANCfc>/ 

CALL  PLTPOF 

XHAP*3.0 

YFAP*?.0 

XL»7.P 

YL=*.* 

M«NCVri 
XS*1 .OF* 

...........................  GRAPH  PMt  - - ------- 

CALL  FIXSCA  (TIM (1) .N.XL.XS1 .XMIM .XMAX1 ,DXl ) 

CALL  FIXSCA  (U3( l> .E.Vl .YSl.YMlNl.YMAXl  ,DY1) 

CALL  COKStA  (L2U1.N.YI  .  YS  1  .  YM  IM  .  YM  AX  1  ,  D  Y 1 ) 

CALL  PLTSCA  (X0AH.YPAR.XMIM.YMIM.XS1.YS1) 

CALL  PLTAXS  (0X1  »0Y  1  * XM INI  . XMAX 1  .  VMIM  , VMAX  1  •  A  ) 

TFMl*XMAXl*nxi/3.0 

CALL  PLT0T2  ( 1 . 0 . T I  M < l)  . U2 <1 ) . h . 0 > 

IJ2N  x  U?(M 

CALL  PLTSYM  (  .  1 . 3M02  >  ,  0 . 0  ,  TEM  .U?  N  ) 

CALL  PLTDT2  ( 1  .O.TIM < 1 ) ,U3 ( 1 ) .N.O) 
l'3M  ■  U3(M 

CALL  PLTSYM  I  .  1 ,3hCR>,0.0.  TEH  *U3. N  ) 

CALL  LAHEI  A  (0  X 1  .0  Y 1  .  XM IM  .  XM  AX  1  ,  VMM  .  YMAX 1  .  XS  .  1 . 0 ) 

CALL  PLTSCA  (XBAB. YEAR. 0.0, 0.0, 1.0. 1.0) 

CALL  PLTSYM  (  .  1 . SVM2 <1 ) . 0 . 0 .  3.0. -0.6) 

CALL  PLTSYM  (.1«SYM1(]),OO.G.  -1.2.1. A) 


CALL  PLTSYM  (,1.PHL0CATI0N>.0.0, 


3.0. -1.0) 


CALL  PLTSYM  ( . 3. 2Hf > . 0 . 0 .  3.0. -1.1) 

CALL  PI  TSYm  l.3.?H]>.n.0t  6.6.-1.1) 

CALL  PLTSYM  (.1.ETA2U)  .0.0.  A.l.-l, 

FbCOPF  (8.50,STEP)ETAn2 
CALL  PLTSYM  (.l.STFF.O.O.  A.fl.-l.l) 

FMCOOF  (P.RO.STEP)ON 

CALL  PLTSYM  ( .  1  , STE P . 0 . 0 .  S.7.-1.1) 


CALL  PLTSYM  ( . 1 . STE P . 0 . 0 .  S.7.-1.1) 

GRAPH  TMO  —  ———  —  —  < 

YH  AM* 12.0 

CALL  PLTSCA (XHAE.YPAP.O. 0.0, 0,1.0, 1 .0) 

CALL  PLTSYM ( . 1 , SY  *(11,00.0.  -1.0, l.«) 

CALL  PLTSYM*. 1.SYM2U)  .0,0.  3.0. -0.0) 

CALL  FIXSCA  ( C  IM  1 )  ,N , Yl . YS* . Y* I M* . VMA X* ,0Y* ) 

CALL  CONSCA  (STC ( 1 ) .N.YL . YS*,YMIN*.YMAXA .DY*) 

CALL  CONSCA  (  TNP  ( 1 )  ,N  ,  YL  «  YS*  «  YM I S*  .  YM  A  X*  ,  l)Y*  ) 

CALL  CONSC.A  (T  AMPL  T  ( 1 )  ,N,YL«  YS*.  YMIN*.  YMAX*,DY*) 

CALL  PLTSCA  ( XOAH , YP AH . XM IN 1 , YM IN* , XS 1 . YS* ) 

CALL  PLTAXS  <CX1,DY*,X*IM,XMAX1,YMIN*,YMAX*.*> 

CALL  LAPELA  *0 X 1 .OY* . XM IM ,XMA X 1 . YM IE* . YMAX* .XS . 1 . 0 ) 
CALL  PLTl'T2(  1 ,0.T1M(  1 )  «C  IM  1  )  «N«0) 

CALL  PLT0T2 ( 1 . O.TIM { 1 ) .STC(l).M.O) 

CALL  PLTDT?  (1.0,TI*'<1)  .  0  AMPL  T(l).N.O) 

CALL  PLTDT2* 1 .O.TIM ( 1 ) ,TNP ( 1 ) .N.O) 

- - .......  STRAIN  PLOTS - - — - 

DC  100  1*1.0  STUN 
J* l *MAXC* ( T- 1 ) 

'VI'  TPA  f  i  ■  |T>.FCCC(  .rliCTJ  ,  *  >.rr  T  ) 

100  CONTIMIF 


MATE 

MAIN 

MAIN 

*>AIn 

GEAPH 
GM  APH 
HE  A  Mb 
GEAHH 
GEAHH 
GPAPM 
HE#  EH 
GEAPH 

geaph 
GEAPH 
GEAPH 
GEAPH 
GEAPH 
GEAPH 
GEAPH 
GEAPH 
GEAPH 
GPAPH 
GEAEH 
GPAPH 
GPAPH 
GEAPH 
GEAPH 
GEAPH 
GEAPH 
GEAPH 
GPAPH 
GEAPH 
GEAPH 
GPAPH 
GPAPH 
GEAHH 
GPAPH 
GEAPH 
GPAPH 
GPAPH 
GPAPH 
GPAPH 
GPAPH 
GHAPH 
GPAHH 
GE  AEH 
GEAPH 
GPAPH 
GEAPH 
GPAPH 
GEAHH 
GEAHH 
GEAHH 
GEAPH 
GPAHH 
GE  A  PH 
GPAPH 
GEAPH 
GPAPH 
GEAEH 
GEAPH 
'■.S.XP- 
GEAPH 


101 


PLOTTING  P(M.  9/*p/75 
IFfNNPF  .FG.OIGOTO  3« 
yp»p  *  4o.o 
00  33  T*).mmHF 
TPAfi  *  YhaP  ♦  10. 0 
IFIVPAP  ,GT.  3  0 . 0 1 C  ALL  PL  TPftF 
IFITPaP  .C'T.  3n.0|VHAO  a  2.0 
FNfOOF  (7,30.«TFP)NPF(I) 

30  FORMAT ( 2HP ( . I3.?H) >) 

CALL  PLTSCA(X8aP.YRAB, 0.0. 0.0. 1.0.1 .0) 

CALL  PLTSYM(.l.$YM?(l) ,0.0,  3.0.-0.F) 

CALL  PLTSYH(.l. STEP. 90.0.  -1.2.2.5) 

CALL  FIXSCA (PD ( 1 , I ) ,N.YL  «YS1 .YPIM ,  YMAX 1 *DY 1 ) 

CALL  PLTSCA  (XHAP.YPAP.XMm  .YMJM  ,XSI  .YS1) 

CALL  PLTAXS(OXI.OYl.XPINl.XPAXl.YPINl.YMAXl.A) 

CALL  LAPELA(CX1.0Y1.XPIM.XPAXI»YPIN1,yi.aX1.XS.)  .0) 

CALL  PL  TOT? ( 1 «  0. T IP ( 1 ) <RO(l,I) , N • 0 ) 

33  CONTINUE 
3N  CONTIMIF 
CALL  PLTPGE 
PFTUPN 

SO  FORMAT  ( F  7 .3, 1 H> ) 

FNO 

Sl'RBOUT  I NF  STRAIN  (X.Y.7,  ETA?,  PM.NETA.J.M 

DIMENSION  X(1 ) . V l l ) ,Z(1) ,  SYM2 (31 ,SYM3 (?) .SYR4 (?) .SYM5 (21 « 

1SVN4<2)  ,>1  (?)  ,X?(?) 

LFVFL?. X.Y.7 

OATA  I/O/. YS/100./.XS/1 .OEh/ 

CATA(SYN2(K) ,K>1.3)/10HFTA?  *  ,10H  (  ,SH  )>/,’ 

?  <SY«3(K) /iontINF  (PICP.10HO«ECOPDS) >/. 

3  (SYM4(K),K«1,2)/1onsTRAIN  <*).|h>/. 

A  <SYM«j(K) .Kal .2) /10HFTA1  COpPO, 5HNENT >/ , 

*>  (SYMA(K) .8»1.2)/inHFTA?  COPPO. 5HNENT>/ 

IF(I.FO.O)GOTO  10 
I»I*  I 

ypabn«ypapn*io.o 

IF  (I.LF.3IG0TC  25 
CALL  PLTPGF 
GOTO  20 

10  CALL  PLTPGF 
X  L  *7 . 9 
YLsN.A 

CALL  r  JXsC  A  ;  X  i  )  )  ..'•«XL.XSI.AyiN,Xv4X,L‘X, 

?0  YFAPNxO.n 
I»1 

25  YPApa?.0.YOAHN 
XP  Apa? . p 

CALL  PLTSCA  ( XHAP.YhAP, n.0.0 .0  a  1 .0, 1  .0 ) 

CALL  PLTSY-  ( . 1 .SYP* ( 1 ) , 90 . 0 .  -1.2. 2. 5) 

CALL  PLTSYM  (  .  1  , S YP 3 ( 1 ) , 0 . 0 ,  3.0.-0.M 

CALL  PLTSYM  (  .  1  .9HL0CA  T  1 0N>  ,  0.0*  3.0, -1.0) 

CALL  PLTSYM  (.3,2Ht>,0.0,  3.P.-1.1) 

CALL  PLTSYM  (.3,?H)>.0.0,  ft. A. -1.1) 

IF ( NET  A .NE • 0 ) GOTO  30 

CALL  PLTSYM  (  .  1  ,NHCLTF.P>,0,0»  A.P.-l.O) 

GOTO  35 

30  CALL  PLTSYM  ( . I ,6HINNFB>,0.0.  N.P.-l.O) 

CALL  PLTSYM  ( . 1 , SYP? ( 1 ) , 0. 0.  A.l.-l.l) 

FNCOOF (H .50 . STEP ) E  TA? 

CALL  PLTSYM  ( . 1 .STFP.O.O.  4.P.-1.1) 

ENCODE (P.50, StFP)pn 1 

LALl  PLTSYM  (  .  i ,  i  T  r  P  ,  0  « ti  ,  *.>*/, -i.l) 
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n3 

GRAPH 

44 

GRAPH 

45 

OF  ALh 

44 

OP  A  PH 

4  . 

GPAPH 

4  8 

GFAPH 

4  V 

C-FAPH 

70 

GFAPH 

71 

GRAPH 

72 

graph 

73 

graph 

74 

GRAPH 

75 

GRAPH 

74 

GRAPH 

77 

bRAPH 

78 

GRAPH 

7» 

GRAPH 

80 

GRAPH 

81 

gpaph 

P? 

GRAPH 

83 

Gf  AHH 

84 

GRAPH 

85 

sthain 

? 

STRAIN 

3 

STRAIN 

4 

STRAIN 

5 

STRAIN 

4 

STRAIN 

7 

strain 

n 

STRAIN 

9 

strain 

10 

stpain 

11 

strain 

1? 

STRAIN 

13 

strain 

14 

strain 

15 

strain 

IN 

STRAIN 

17 

STRAIN 

IP 

STRAIN 

19 

STRAIN 

?P 

strain 

?1 

STRAIN 

STRAIN 

?3 

STRAIN 

2* 

STRAIN 

?5 

STRAIN 

24 

SThAIN 

27 

STRAIN 

28 

strain 

29 

strain 

30 

strain 

31 

strain 

32 

strain 

33 

STRAIN 

34 

STRAIN 

35 

STRAIN 

34 

strain 

37 

STPAIN 

38 

STRAIN 

35 

STRAIN 

40 

s 7  R  h  ii- 

-  i 

xi  <n»?.o* 

XI  (2)  ®3.S70 
X? (?) *-l .3 

CALL  PLTOTS  (l .O.X] ( 1) .X2(l) .2.0) 

CALL  PLTSY*  (.l.SYMS(l) .0,0.  3. 7, -1.3) 

X2 ( l ) *-1 .5 
X2(2)=-1.S 

CALL  PLTPTS  (A  .1)  .X 1  ( 1 )  .X2  <  1 )  .2.0 ) 

CALL  PLTSYM  (.  1  .SYMf.  (1 )  .0.0.  3. 7.-1.?) 

CALL  F I XSC A  <Y(J> .N.YL.YSl.YNIN.YMAX.OY) 

CALL  CONSCA  (7  <  J)  .N. YL.YS1 . YMIN.YMAX.DY) 

CALL  PLTSCA  (XPAR.YPAP.XMlN.YMIN.XSl.YSl) 

CALL  PLTAXS  (CX  .OY . XM IX. .  XMAX  . YM  IN .  YMAX  •  A ) 

CALL  PLTDT2  < 1 . 0  .  X ( 1 ) . Y ( J) .N, 0 ) 

CALL  PLTDT?  1 A  .  O.X  ( 1 ) . 7 ( J) .N. 0 ) 

CALL  LABELA  < DX . OY . X M I X . XM AX . VM IN , YM AX . X S * YS ) 

RETURN 

*0  FORMAT  (F7.3.1HX 
FNO 

SUHPOUTINF  FIXSCA  (X.NPTS.SI7F.XSCALE.XMIN.XMAX.-0X) 
DIMENSION  X(l) «T (2) 

MUST  BFAO  IN  P I X SC A  TO  PUT  X  IN  LFVEL  ?. 

LFVFL2.X 

I)  X-Al  INFAP  ARPAY  OF  NUMBERS  IN  DATA  UMTS 
I)  NPTS-THE  NUMBER  OF  X  VALUES 

I)  THE  LENGTH  OF  X  DIFFUSION  OF  THE  t-RAPM  IN  PLOTTFH  UNITS 
9)  XSCALF-THt  SCALE  IN  DATA  UNI T S/PL OTTER  UNIT 
P)  XMIN-THE  AOJUSTE  0  MINIMUM  IN  OATA  UNITS 
t,)  XMxx-THP  ADJUSTED  MINIMUM  IN  OATA  UNITS 
9)  nX-THF  DELTA  X  FOR  AXIS  IN  OATA  UMTS 
LOGICAL  CONT 

cont*. false, 
goto  100 

FNTPY  CONSCA 
CONTm.TPUE. 

100  TX**!*X  (1  ) 

TXMA*X(1) 

DO  120  I«1,NPTS 

IF  (X(I)  .(iF.TXMI)  GOTO  110 

TXMI«X(I) 

GOTO  120 

110  IF  {X ( I)  .1 E.TXMA1GOTO  120 
TXMAax <  1  ) 

120  CONTINUF 

IF (  .NOT .CUNT ) GOTO  1A0 

IF  (TXH.lT.XMIN.OR.TXMA.GT.XmaXIGUTC  12? 

TXM I «XMJ N* APS (XM IN) #0.000001 
TXMAsXMAX-AHS (XMAX) #0.000001 
GOTO  1  A  0 

12?  IFITXMl.LT. XMIMGDTO  130 

TXMI«XMIN*AHS (XMIN) *0.000001 
130  IF ( TXMA.GT . XMAX ) GOTO  1*0 

TXMAsXMAX-APS (XM AX) *0.00 0001 
1*0  OIFF.TXMA-TXMI 

IF  (DIFF.LT. 0.0 00001  )  01 FF*0. 000001 
FNCOPF  (20.  ISO.  T)  DIFF 
ISO  FOPMATdPF  12. S) 

OF  COPE  (20.  100.  T 1  COEF.  IEXP 
too  F OPMaT (FP.5.1X.F3) 

170  IF  (COFF.GT.2.01GOTC  100 

rsi.  T4.fi .  i 

r-nro  200 
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IflO  IF  (COFF.c.T  .4 . 0 1  GOTO  too 
f!f  LTA«0»? 

GOTO  200 
ton  CEL  TA*n,o 

200  HALT A*r>EI  T A*  1  0 «0** IF XP 
OX*OFl TA 

XK!f.«AlMT(4«S<TXKI)/Cei  TA)*l)FLTA 
IF  (TXKI.LT.il >  X*TF«-  <X*IN*OfclT A) 
XKAXaAIKiT  (  AHS  (TXKA)  /rFLTA)*Ct'LTA 
IF  (TXMA.t T . 0 ) GOTO  210 
XMAX*XMAX ♦DFl  T  A 
GOTO  220 
210  XKAXa.XMX 
220  XSCALF«(XKAX-XKIn) /SIZt 
HETHHN 
FNO 


FIXSCA 
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APPENDIX  E 


EXAMPLES 


Two  examples  are  given  in  this  appendix.  These  were  originally  jobs 
run  on  the  BRLESC  computer,  now  defunct,  and  transferred  to  the  CDC  CYBER  70 
Model  76  computer  to  check  the  program  there.  The  two  jobs  test,  and  hence 
illustrate,  many  of  the  available  options  and  internal  variations  of  RPSL1D. 
For  both  examples,  the  complete  SCOPE  2.1  batch  job  is  listed  and  a  sample 
of  the  output  is  given.  The  sections  of  the  batch  job  are  separated  in  the 
listings  by  question  marks  which  listed  for  the  multipunch  end  of  record  and 
end  of  information  signals. 

E.l  Example  1:  Pressure  Loaded  Flat  Plate 

The  first  example  involves  finding  the  response  along  the  shorter 
symmetry  line  of  a  flat  4  by  15  inch  soft  aluminum  plate.  (Results  from 
the  standard  REPSIL  and  the  slab  symmetric  version  of  RPSL1D  are  very  similar 
for  this  elongated  plate.)  This  example  demonstrates  the  introduction  of  a 
new  PRESS  subroutine,  the  use  of  the  optional  coding  PLOTP,  and  the  activating 
of  Gaussian  integration  and  damping.  Both  the  new  PRESS  routine  and  the 
option  PLOTP  are  inserted  by  making  changes  in  the  program  through  UPDATE; 
both  also  require  input  data  cards.  With  the  catalogued  INGEOM,  Gaussian 
integration  is  automatically  used  with  slab  symmetry  if  IB  =  0  on  input  card 
14.  Damping  was  activated  from  time  At  on  by  setting  MDAMP  =  0  on  input 
card  5.  The  damping  option  was  used  here  to  help  simulate  displacement  of 
the  plate  by  quasi-static  pressure  loading.  This  is  not  the  usual  use  of 
damping.  Damping  is  a  device  normally  used  to  bring  a  responding  surface 
to  rest  at  its  final  equilibrium  configuration  after  the  loading  is  removed. 
Since  subroutine  DAMP  is  not  entered  until  cycle  one,  when  NCYCLE  =  1,  the 
subroutine  did  not  shut  off  pressure  loading  as  it  will  if  entered  with 
NCYCLE  =  MDAMP.  This  run  terminated  at  cycle  3546  with  a  displacement  of 
. 66  inches . 

A  listing  of  the  batch  job  deck  and  samples  of  output  for  Example  1  fol¬ 
low.  The  first  part  of  the  job  deck  is  the  SCOPE  2.1  control  statements  to 
set  up  and  run  RPSL1D  and  the  plotting.  The  next  set  of  cards  is  UPDATE 
directives  and  new  FORTRAN  statements  for  RPSL1D.  The  input  data  cards  for 
RPSL1D  are  next.  (The  input  and  output  for  this  example  are  in  the  now 
forbidden  pound -inch -second  system.  If  this  is  discomforting,  pretend  the 
units  are  SI  units  for  some  very  exotic  material.)  After  the  data  there  are 
UPDATE  changes  for  the  plotting  program,  including  changes  for  the  PLOTP 
option  and  a  change  to  permit  up  to  6000  time  cycles.  Finally,  there  is  a 
data  card  to  control  plotting. 

Following  the  job  listing  are  sample  listings  of  output.  This  starts 
with  a  summary  of  the  input  and  the  initial  Cartesian  coordinates  as 
prescribed  by  INGEOM.  All  the  output  should  be  inspected  to  uncover  errors, 
but  these  two  sets  should  receive  very  close  attention. 
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The  output  at  cycle  3000  for  this  example  is  listed  next.  This  includes: 

•  Displacement  increments,  coordinates,  and  pressure  at  each 
mesh  point  including  the  virtual  external  point  at  N=1  and 
the  symmetry  point  at  N=13. 

•  The  LMAT  matrix.  A  "1"  signifies  plasticity,  and  a  ”0"  denotes 
an  elastic  response,  on  the  current  cycle.  For  each  column 
labeled  N,  there  are  two  rows  for  each  integration  station  K. 

The  first  row  is  for  mesh  point  N.  The  second  row  is  for 
midmesh  N. 

•  The  components  of  the  surface  normal  at  each  mesh  point,  the 
2  columns  on  the  left,  and  at  each  midmesh. 

•  Surface  strains  at  all  mesh  points.  For  each  mesh  point, 
the  strains  on  the  upper  surface  are  listed  first,  then  the 
strains  on  the  lower  surface  are  listed.  EC1  and  EC2  are 
the  elongation  strains  in  the  S1  and  E2  directions,  respec¬ 
tively.  Ell  and  E22  are  internal  covariant  components  of 
strain.  For  slab  symmetry,  EC1  and  Ell  are  zero. 

•  The  tabulated  output  for  the  NSTRN  prescribe  surface  strains 
accompanied  by  prints  of  extreme  strains  and  maximum  deflec¬ 
tion  for  the  current  cycle,  and  for  the  entire  run. 

•  The  energy  balance  summary. 

In  our  listing,  this  is  followed  by  three  prints  noting  the  use  of  the 
kinetic  energy  annihilating  procedure  (a  part  of  damping  which  stops  motion 
at  kinetic  energy  maxima)  and  finally  a  print  of  the  prescribed  surface 
strains,  etc.  at  cycle  3500. 

The  occurrence  and  frequency  of  this  output  is  controlled  by  input  cards 
8,  9,  and  10. 

Figures  E.1.1  through  E.1.4  are  samples  of  plotting  output.  For  most 
publication  purposes,  it  is  desireable  to  do  some  relabeling,  but  these  plots 
have  not  been  retouched.  Figure  E.1.1  shows  displacement  at  selected  time 
steps.  Figure  E.1.2  shows  the  displacement  of  the  selected  central  mesh  point 
and  the  energy  balance  as  functions  of  time.  The  scallops  in  these  plots 
occur  when  the  stopping  procedure  was  exercised.  This  energy  balance  plot 
shows  the  four  possible  curves.  These  are,  from  low  to  high:  kinetic  energy, 
kinetic  energy  plus  elastic  strain  energy,  kinetic  energy  plus  elastic  strain 
energy  plus  the  work  in  damping,  and  the  total  energy.  Figure  E.1.3  is  the 
first  two  of  the  six  strain  plots  requested  for  this  run.  Figure  E.1.4  is 
a  plot  of  force  per  initial  unit  area  at  mesh  point  3.  The  pressure  increased 
linearly  to  a  maximum  of  80  psi  and  then  remained  constant,  but  the  area  in¬ 
creased  throughout  the  entire  run. 
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Listing  of  the  Batch  Job  Cards  for  Example  1 


Wf'f  T>'«STM*7.TM'.  MPl.KUF  BLN  If  1/3/77.  TNST  FAUSS.  PLC  TP •  LINPC*. 
»n:i  i*  T. ****«.  *Of  T!»  4N  •*,’('B»  A3P7<5 

attacp <oL r t|  .cFSL1l'.To=.Ji  » ) 
upr  atf.p  . 

FTMfl.I.SI  *0.P*0) 
vao.OFF. 

L«r. 

FKTTO.M  OATA  TO  PLOT  T  APf 

PFu  TNO (PLOT) 

CFTliFN .COMPJl  £ . 

qfthfn.loo. 

PFTlIPN  ,01  OPL  . 

ATTAO  IOLOPL  .PPSI  1 1'PL  OT «  Tt'sj('W) 

pf fin. attaop.pl ctlTp. 

IIPOATF  (F) 

FT*  (A,  I  .Fl  aF.SI.  *0  ,K=n  ,LfP*I ) 

PAP. OFF. 

LOO. 

F  *  T  T  ( t  ) 

“FO U .PLOT.CAl COMP.TAPF 1 7. 

7 

ofnFA.T  PI. OTP 


*  T 

PPSL 10.F 

C 

PI  OTTIAO  PIN).  <5/7n/7p 

PL  OTP 

»T 

»«!*.?» 

OOPFON  NNPF.NPFI  Pl.enATI  0) 

PLOTP 

°  T 

*=takT.?4 

PFAPI*. 1 1  0>NNPP  .  I>.FE(I)  •  1*1  .NNPF) 

PLOTP 

OT 

Fr.ATA.).7 

PPTTF (f  PLOT1M.PF . INPF ( I)  . 1*1 .NNPF ) 
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*T  . 

POATA.IO 
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O  T 
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1  •  (PI'AT  1 J)  .„•!  .NNPF  ) 

PLOTP 

*  T 
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Co  AN  T  =  1  .>'NFN 

PLOTP 

w  =  l>  Pc  (  i  ) 

r  L  i>T  P 

/.c  PDATITlxPIJ) 

PLOTP 

*1 

C|'4T#,,F 

1  ,  (PC  AT (J)  ,w«1  ,M»BF ) 
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o  Tr 
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or 

BCF  C<=  .4  ,P|-P  «!«  .  1  4 
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L INpPS 
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on  MP1PI.T.40 

1  .  (PPftT  ) 

»tr  Penn? 

or-  "A  Jp|.7.**o  [p  .q 

1  T  T  p  (  \.~  (  pm?),  till  Pl  o? )  •  C1M  P  0  0  r  )  t 

?  Tvt  (  pn? )  , iiAfPl  t  <  p(.o?)  «FPSS1  OPol?)  »€PSS?  OPO 1?  1 

C<lf»MiK/M  0TP/PO  (  POOP.O) 
oo  pPJPLT.in 


STO (  POO?). 


*003  II  $44(1  I'M 


5  c  • 


IS  *5 


S  £& 


s  »i 

i  J: 

9  ft  W 

;  ^ 

•  Jc 

ftft 


•  (hi 

4 V  5 
ft  ••* 

t  :« 

*  V 

•  wi 


:r  •• 

ft  *  ft 

■s  *  i! 

*  ►  ft  * 


i  re? 

«  «<u 


t  ::  £5"‘ 

r*ss  Is 

yT  |  | 

ft-ft  V 

«  ►A** 

»i£  ii-s 

srss  g;rr 

—  *s-* 
:tf 
.  •  t** 


*  w  —  ft  ft 

c  *- 

e  a  —  *  ^  ft  ft 

Mi  s  t 

ft  *  If  _  8  •*  “ 


s:r«t  •  : 

e  «K.fc  -  I  t 


2§*"*'*  •»  —  w 

t.ssi  s  s  5 

s-:::  i**,s 

—  •  •*««»> 


t  ft  ft  * 
hft  •  ft 


••ft  «  •» 

ft»  ft  M 

•  ft*  . 

I-  :  I 

—  ft  »> 

Man  ft  #•!»• 

is  -  t?t:  i 

“i  $;£us  : 

•  *•••••  * 

«►  <•••••  • 

«  —*••••  ft 

ft  «*!»•*•  —  • 

—  —  ft  ft  m  ft  »  ft  ft« 

ft  K  ►  •  t  •  •  |  < 

•  C  ft*  ft  ft.  —  ft  ft  4 


«•  ft  ft  ft  ft  ft  ft  ft  ft 

•  •  «■  III  ft  •  •  •  • 

ftftftftft  ft*  ft  W  ft  ft  ft 

X  ft  ft  •  o  — —  ft  »  ft  ft  ft  * 
ftft  ft.  ft  ft  ft  *)•••• 

_  —  _  ^  SI  —  ft  ft  ft  ft 
ftftftftO  ftft  ft  ft  ft  ft  ft 
«  •  5  ft  ft  ft«r.  —  •  r  ft 

•  *  •  ft  •  v«ft ft  ftft 

r.ft  r.  *  ft  «  *>, 
*  _  ■*  X  ftft  •••• 

•■ft  ft  —  ftft  ft  ft  —  — 

ftftMft  ui 

5  ft  —  ft*  ft 

ft  —  ft  —  ft  ft 

^ssSftS 

«— w£©5  ^  •.**** 

sssscs 

rr^ 

ft  ft  pi 


Selected  Output  from  Example 


IHMM  prikl  MINI  |  INI  WIN! 

.,iimmhmjij«im*m  -.iiMMimiinmiiii  i. 


iWfW«wiM«wiw*.*<vfe 


xirr>r>**c«v«*> 

•  if 

•  *»»-*•  -0  x  e  ♦ 


•*•*■:*•  ir  «v  —  *  •  » 


*•  —  **-«»  — 
»  f  «  <  •  * .  <r  k  c  u 
9  —  —  —  —  N  IV 


#*  —  If.  *»••«*  —  IT  — 

r.  £££••*•£  —  £ 

iO  iv*  —  —  «w  if.  *  *► 


NC«CNCC4C« 

o-  o  *  •  »  if)  **>  •  m  ♦ 
*  iv  •  •  —  *  *  —  rj  *• 
fvive  —  n^AfwriK 
iv  r*.  •  «  ir.  ir  c  f  tf  < 


•  • 


i  •  • 


•  I  I  I  •  II  II  I  I  I  — 


ir»  ir  *>  rift  «r»  #i  ir. 

Bccceeec 
i  I  I  (  •  •  »  i 


*— r*«n  —  «iv« 


a 

*» 

— 

•f. 

r 


l  l  i  l  •  •  l  i  •« 


•  •  i  i  i  i  •  «  i  i  i 


i r«*»«<vi»*»e0» 
lT0^4«^C«NCN 
c  ir.  if.  ^  •  —  ♦  c  c  r  & 

e  r .  tr>  ^  <  rv  —  —  — 

**■  <C  iv  rv  •  *  *  iv  ir  fv 

c  *  —  —  c  • 

IV  »  C  »  *  it  I  K  i  c  1 

«  ^•lw»«*»fv«» 
«  r.  c  rw  fk  •  r*.  —  ir  e  / 
(v«e^*»?— oe. 
ifetfcifr.  iififtff 
creec*rv*-c— c 
t  — ^r.  —  —  — 

^  i  i  i  i  i  i  i  i  i  i  • 


t*NP«lTil»C0fNr, 

—  -  -  —  X 


111 


www  «ii»in«..  ww«rt«wi>.  itti.ici>niiim 


cnir(«« 


Z«««5*2m«  22222E2CJJSSS 

»  C  •  0  —  «  •  9  p-  n 

i.‘  kn»  «Nt.» 

«•>»£  —  —  •»  #■  9  « 

*>  iw  *v  *  «r.  t  n«iw« 
<  <  If-  «  «  «  •  9  *  c  —  9  •  < 

n  n  ^  •  t  «»NNC»f  »•  «r.  at  *  »r.  *>«nn 

9  9  •*  pi  • 

•  *•••••••••••••••••••• 


c 

* 

s 


I 

t 


« 

9 

*■ 


e  e  r  eec 

kizttz 


2  «  9  «w  s  S 

«  •  9  «  *1  IT 
m  tfi  *> 


! 

E 


« r> 

* 


«  M  M  (W  «  M 

•  I  I  I  »  • 


■•fwMr  iDjn 

4>  p»  P»  9k  P>  « 
•  •••«• 


•22*25 
•  •••»• 


as 

s 

• 

•  «•  —  IW  —  •■ 

«  1  1  1  1  • 

m 

F. 

■ 

« 

9 

r>  —  *>  ir.  *  F 

Ito 

a 

IW««N9  i 
«  <vn  «  »  pi 

W 

IT. 

* 

?s?sss 

9 

U 

X 

« 

t«Sr?5 

I  •  •  • 


»•««••«• 


•  •  •••••• 

•  «w  •  9  «  ri^r  kr  r.r 

•  <  9  •  m 


•.:tv.v.n 

^  ■  t  •  ~&S22S2SS 

^S2S2j£.  222 

2  222222222 

«  »«••••• 

2  222222222 
*  SSSaSSgtf* 


222222222 


to  »  ^  rv  c  —  rv 


fw  *  •  r  «  *  c  t 


9 

& 


IIDIfWI  I.M.l 


CYCLE  O 

CeFL£CT?C$  l  -COO 


SCPLE  1/  O 


Figure  E.1.1 
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Plots  of  the  reference  curve  at 
selected  times. 
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Figure  E.1.3.  Strain  histories  at  the  edge  on  the 
top  and  bottom  surfaces 
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.4.  Force  per  initial  unit  area  at 
mesh  point  B. 
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This  example  uses  the  beam  option  of  RPSL1D  and  utilizes  all  seven 
of  the  options  listed  in  Section  5.  This  program  was  designed  to  illustrate 
that  the  test  strain-time  records  at  two  locations  along  a  steel  rod  pene¬ 
trating  a  steel  plate  could  be  reproduced  by  varying  the  force  and  the  rate 
of  erosion  on  the  end.  Two  quite  different  tables  of  erosion  and  force  were 
used,  and  both  returned  strains  that  were  close  to  the  test  values  at  three 
locations.  This  example  used  a  fast  erosion  rate.  Mesh  points  2  and  3  were 
eroded  by  time  cycle  5. 

The  batch  job  cards  and  some  sample  output  are  listed.  The  batch  job 
for  this  example  is  rather  long.  The  first  section  of  the  run  stream  is  the 
SCOPE  2.1  control  statements.  The  next  section  is  the  UPDATE  directives  and 
new  FORTRAN  statements  for  RPSL1D.  Breaking  this  section  into  the  nine 
groups  with  different  UPDATE  card  identifiers  we  have:  PLOTP,  SHR3/1,  MSQSVS, 
BSTRS,  APLFRC  (including  a  subroutine  ENDFRC  which  is  adjusted  for  erosion), 
EAPFRC ,  11/2/76  (a  modification  for  EAPFRC) ,  INVELC  (an  insert  into  INVEL 
which  creates  a  uniform  initial  velocity  toward  end  1),  and  ERODE. 

The  next  section  in  the  job  cards  is  the  input  data.  The  input,  and 
output,  for  this  example  are  in  SI  units:  meters,  seconds,  kilograms, 

Pascals,  etc.  This  leads  to  some  awkward  scaling  on  both  input  and  output, 
but  assures  a  consistent,  unambiguous  system.  Along  with  the  usual  input 
data,  this  section  contains  a  long  table  of  strain  vs.  stress  for  the 
BSTRS  option,  input  to  approximate  a  solid  cylindrical  beam  for  INGEOM,  no 
input  for  PRESS,  a  card  with  initial  velocity  for  INVEL,  three  tables  for 
ENDFRC,  and  a  table  for  ERODE. 

The  next  section  is  UPDATE  changes  for  the  file  RPSL1DPL0T  for  options 
MSQSVS  and  PLOTP.  Finally,  there  is  an  input  data  card  to  control  the 
plotting . 

The  tabular  output  listed  is  two  of  the  seven  pages  which  reflect  input 
(most  of  the  omitted  tabulations  are  simply  copies  of  input  tables) ,  and  the 
output  at  cycle  130.  The  output  at  cycle  130  includes: 

•  The  displacement  and  coordinates.  R  and  DR  are  effectively 
zero. 

•  The  output  from  SHR3/1.  The  moment,  MS,  and  shear,  VS,  are 
zero.  Only  axial  force,  QS,  is  significant. 

•  The  table  of  plastic  stress  occurrence. 

•  The  components  of  the  unit  normal  vector.  The  rod  has  not  bent. 

•  The  surface  strains. 

•  The  energy  balance. 


Six  of  the  output  plots  are  shown  next.  (An  isometric  plot  of  the  rod 
and  plots  of  MS,  QS,  and  VS  were  made  at  cycles  70  and  130.  The  isometric 
plot  is  uninteresting,  and  MS  and  VS  are  zero.)  A  plot  of  QS  vs.  mesh  point 
number  at  cycle  130  is  included  to  illustrate  output  of  the  MSQSVS  option. 

We  also  include  plots  of  several  functions  of  time:  the  displacement  of  mesh 
point  S,  strain  plots  at  19.26  mm  and  39.67  mm,  a  plot  of  force  on  end  1,  EFZ1, 
which  was  stored  in  P(53),  and  the  erosion  at  end  1  which  was  stored  in  P(54). 


Listing  of  the  batch  job  cards  for  Example  2 


WOPTM.STMF/.T30.  OUPLICATF  PUN  OF  6/A/77  AND  10/12/77.  ERODING  PCD. 

ACCOUNT .  W0«TMAN,P3n9,X3979 

ATTACH, OL  DPL  .PPSLIO.lDaJOW. 

UPOATE.F. 

ETNIA.I.SLaO.RaO) 

MAP. OFF  . 

LOO. 

EXIT  t U )  DATA  TO  PLOT  TAPE 

REMIND (BLOT) 

PFTURN. COMPILE  . 

PF  TOWN , LOO . 

RFTUWN.OLOPl . 

ATTACHIOLrPL.HPSLlOPLOT. I0sjr») 

MFC’ IN.  ATTACH,  PLOTL  IP. 

UPDATE (F) 

FTN(A,I.EL*F.SL=0.Ps01 
MAP. OFF. 

LOO. 

EXIT  (111 

UFO I N .PLOT , C  ALCOMP . T  APE  13. 

*  I l'i ENT  PLOTP 
*1  PPSL10.S 

C  PLOTTING  PIN).  9/30/7S  PLOTP 

•  I  *<AIN.?e 

COMMON  NNPE.NPEl  9),P0AT<  9)  PLOTP 

•I  START. ?A 

PFA0(9.110)NNPF. (NPE ( I  1 , lal.NNPF)  PLvTP 


*1 

P0ATA.13 

WPITF (NPLOT)NNPE. ( NPE (I 1 , I» 1 .NNPE ) 

PLOTP 

*1 

PCATA.IO 

DO  ?H  lal.NNPF 

PLOTP 

JaNPE(I) 

PLOTP 

?G  POAT ( I 1 »P ( J) 

PLOTP 

•  I 

PDATA.R7 

1  . (POaT(J) .usl,NNPF) 

PLOTP 

•  I 

PPATA.42 

00  AH  1*1 .NNPF 

PLOTP 

JaNPF ( I ) 

PLOTP 

an  POAT(IlsP(J) 

PLOTP 

*1 

PPATA.4S 

1  . (PrAT ( J) , J»1 .NNPF ) 

PLOTP 

*  I OF  NT  SMP3/1 

*1 

PPSLIO.S 

C 

FSTIMATF  SHEAR  FOWCE,  VS.  ANO  MS  ANO  OS. 

3/1/7H  (FOW  A  REAM) 

SHP 

3/1 

*1 

MAIN.RM 

COMMON  AMS ( 103) .MS ( 103) »VS( 103) 

SHW 

3/1 

*1 

RESULT. 4H 

c 

MS  AT  FPFF  END  IS  ?EW<J  U/V’M.  CHANGE 

V I TH  APLFPC. 

SHP 

3/1 

AMS(N?P)  a  0.0 

SHP 

3/1 

•  I 

RESULT. Sa 

AMS(N)  a  A??*E  ?2 

SHP 

3/1 

•  9 

WFSUl  T.H? 

<S  •••'•!  *  '  i-S  ,  «,.-.)  -  l»SfN)  )  /  iCF?  J?*eWA' 

SE  S 

3/1 

POO  (1S<M  *  A2?« (C?2-P«2?»F22) 

SHP 

3/1 
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VS. MS. OS  PM  1 NT  CONTROLLED  BY  COOIN6.  SHK  3/1/76. 

TO  m P I N T  EVERY  10  CYCLES. 

IF  [MOD  (NCvCLe  .10)  .t'O.  0)  GOTO  Ft? 

!F(»'CvriF  ,eo.  70  .09.  NCYCLt  .KO.  1301G0T0  ?1? 

GOTO  21B 

?10  FOPMflT(//t  CYCLE  N*.«-*.,MS*.13Y.»VS».13X.»«S».13F.,A22,.12X. 

1  •*??• .1?X. *L?2* . 1?X, *R2?» ./) 

?11  FORMAT (TM.t3.lP7E 15.7) 

21?  I F  !  K  .FO.  NIP  (WRITE (6.2101 

WRITF<6.21 l)NCYCLE.N.AMS(N) .VS  IN) .OS(N) • A?? . BM22 .G?2 .F?2 
OUTPUT  OK  MSIN2B)  */5/76 
I F ( N  .EO.  N2P-l)wPITE(6.211>NCYCLE.N?P,AHS(N2R> 

?1P  CONTINUE 
IDFNT  MSOSVS 
I  RPSLIP.S 

PLOTTING  MS.  OS.  AND  VS.  3/19/76 
I  P0ATA.S2 

WRITE  INPLOT)  (AMS  IN)  .OS(N)  .VS  IN)  .NsNlH,N2B) 

IOENT  PSTRS 
I  PPSL1D.5 

NEW  UNIAXIAL  STRESS-STRAIN  CURVE  FOR  BEAMS.  A/1/77 
NOT  COMPATIBLE  «ITH  STREN  OPTION. 

0  MAIN. 5 

3  S 161  (103*6) .61 G? (103.6) .LMAT  (103.6) * 

*0  “AIN.?* 

COMMON  SNRM (103) * SNKM (103).SIGlM(103.f) .SIG2MI 1 03 . 6 ) .LM ATM ( 1 0 3 .6 ) 
*1  MAIN. ?p 

COMMON  tPSP (103.6) .EPSPMI 103.6) .FPLBI 103.6) .EPLBM ( 1 03 .6 )  . 

2  EFR(li)O)  .SSR(IOO)  .NEST.EPS7 
•I  PPSL10.60 

EPSP  (N.K)  *  0.0 
EPSPM(N.K)  a  0.0 
EPLB  (N.K)  a  0.0 
FPL PM ( N , K )  a  0.0 
•O  START. 16. START. ?0 
00  700  I«J.  99 
REAO(S.70?)EFR(I)  .SSP(I) 

IF (FFP ( I )  .LT.  0.0)  GOTO  701 

700  CONTINUE 

701  NFST  a  I 
NSFL  *  1 

SSP(I)  a  SSM(T-l) 

FER(I)  a  1000. 

SSH(l)  a  SIGZ 
EPS7  a  SIGZ/F 
EFM ( 1 )  a  EPSZ 

WPITF (6.703) ( I.EEB ( I) .SSB ( I ) . Ial .NFST) 

I F ( SSM (NFST)  .GT.  2 . 0*S IGZ ) WP ITF ( 6 t 70* ) NFST 
70?  FOPMAT ( ?F 1 S • 7 ) 

703  FORMAT(/3X.*STPESS-STRAIN  TARLE  FOP  PSTRS * /3X . *  I • . 6X . * F PS ( I ) • . 9X . 
?  «SIG( I ) * ./ (IS. 1P?F IS. 7) ) 

70*  FORMAT!//*  «*«WARNING***  STRAIN  FNPHGY  SUSPICIOUS.  SIO(«.I3. 

?  •)  .GT.  ?*S IG7  * // ) 

*0  HMSTRS.73.RMSTRS.30 

C  CHANGFS  IN  PMSTRS  FOR  BSTRS  */l/77 

IF(IP  .FO.  1 ) GOTO  10 
PEPS  a  FRSPM(N.K) 

HEPI  a  EPLPM(N.K) 

GOTO  11 

IP  BEPS  a  FPSP (N.K) 

RFPl  a  FPLH(N.K) 

11  PFP'  M  =  MFDL 

IF(RFPLM  .GF.  0.0)  GOTO  13 


9/15/76 
Tt  ST 
SEP  3/1 
St-R  3/1 
SEE  3/1 
SHE  3/1 
SEP  3/1 
SEP  3/1 
SEP  3/1 
SEP  3/1 
SRP  3/1 


KSCSVS 
K  SOS VS 


tSTRS 

PSTRS 

PSTRS 

PSTPS 

PSTRS 

PSTRS 

PSTRS 

PSTRS 

tSTRS 

PSTRS 

PSTRS 

PSTRS 

FSTMS 

PSTRS 

PSTRS 

PSTRS 

PSTRS 

PSTRS 

PSTRS 

PSTRS 

FSTRS 

PSTRS 

PSTRS 

PSTPS 

PSTRS 

PSTRS 

PSTRS 

PSTRS 

PSTPS 

PSTPS 

PSTRS 

HSTRS 

PSTRS 

PSTRS 

PSTPS 

PSTR-3 

pstrs 
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STC-P2I  *  -S102PI 
0EPSP2  *  -OEPS?? 

reps  =  -reps 

«ERl  =  -PERL 

13  PEPS  a  RE°S  ♦  PEPS ? 2 

IF (PE°S??  .GF.  O.niGPTP  14 

IF ( PFRS  .GF.  RFPL-EPSZIGOTO  IS 

IF ( REPS  .GT.  -FPS71GPTP  16 

CALL  PVOINTI  -HtPS.SIGP2.EFP.SS8.NFST.?) 

SIG?2  =  -SIG22 
GOTO  17 

14  IF (PEPS  .LE.  BEPL *EPSZ I  GOTO  15 

TALL  OVPINT  <8FPS,SIG??.EEH,SSB.NEST.2> 

PEPl.  *  8FPS  -  EPSZ 
L»NK  s  1 
GOTO  1H 

15  SIG22  =  STG?2I  ♦  E*PEPS22 
GOTO  18 

16  CALL  OVOINT <PEPS«2.0«EPSZ.SIG??.EEP.SSP.NEST,?> 

S  IG?2  =  SIG22  -  2 .0*5 IGZ 

17  LMNN  *  1 

8FPL  *  PEPS  «  EPSZ 
IP  IF (HFPLP  .GE.  0.0) GOTO  19 
S  IG?2  =  -SIG?? 

HFPS  a  -PEPS 
PEPL  =  -HEPL 

19  I F ( I P  .FG.  1)  GOTO  20 
FPSHM(N.K)  a  mFPS 
EPLPM<N.K)  a  PERL 
GOTO  21 

?0  EPSB(N.K)  a  PEPS 
FPLh(N.F)  a  PFPL 
21  CONTINUE 
*  I OFNT  aplfpc 
•I  PPSL1D.5 

APPLIED  FORCE  (ANO  MOMENT!  AS  FNO  CONDITION  4.  4/12/76 

APL  FBC  ASSUMES  SMP  3/1  IS  INCLUDED  IN  THE  DECK. 

the  applied  force  change  to  n-D)  repstl  for  beams  reouihfs 
a  user  suppoutine.  enofhc.  to  supply  end  FORCES  AND  MOMENTS. 

I  MAIN.2B 

COMMON  EF91 .EFZ1.EM1 .FFR2.EFZ2.EM2 
I  START. 54 

IF (IPCF1  .FO.  4)N1W*N1P 
h  START.  230 

I F ( I BCF 1  .LT.  5  .AND.  I8CE2  .LT.  5 ) RETURN 
D  GPAD.9.GPAD. 1 0 

I F  (  I  BCF  1  .NE.  4  .OR.  N  .GT.  MB1P0T0  H 

IPCF  1=4  ANO  NsMR.  USE  FORpARD  DIFFFWFNCES 
9?  a  PT0?«(-3.O*  P(N)»4.0*  H<N.]>-  R ( N *2 ) ) 
z?  s  http* (-3.0*  Z(N).4.n*  zcn.u-  zin.pi  > 

DP?  a  PTU2*(-3.0*CR (N) ♦4.0*DR(N»1 )-0R(N*21 ) 

07?  a  PTD?* (-3. 0*C7 <N) .4.0*07 <N»1 1-DZ (N«2) > 

»??  »  RD22®(?.0»  RCNJ-5.0*  R(N»1).4.0*  R(N'»2)-  H(N*31) 

Z2?  a  PD??« (?. 0*  Z (N) -5.0*  7(N« 11.4.0*  Z ( N*2 ) -  Z(N*3)) 

OP??  3  RD??« (?.0*PR (M -5,n*DR ( N . 1  1 »4.0*DR <N*2)-Dfi C N ♦ 3 » ) 

DZ??  a  RD2?« (2.0*CZ (M -5.0*07 (N.l ) «4.0*DZ (N.?)-OZ (N*3) ) 

GOTO  ?0 

P  IF ( IflCE?  .NF.  4  .CP.  N  .LT.  N2H1GOTO  9 

IPCF?a4  AND  NsN?P.  USE  BACKRARP  DIFFERKNCFS. 

GRA0.34.GPA0.35 

IFdMC.Fl  .NE.  *  .CP.  N  .GT.  hlR  (GOTO  18 
C  T  PCF 1  *4  1N0  '.SNIP  <»1/?'.  USE  NOF-LFNTual  t  OOI-.T  5ND  I:  I FF  FRF  NCCS 
P??  =  RP??M»(3.n*  W(N)-7.0*  P<N»l)«S.n*  PIN.?)-  PIN. 3)1 


PSTRS 

PSTRS 

RSTRS 

FSTRS 

hsths 

RSTRS 

HSTRS 

PSTRS 

PSTRS 

RSTRS 

RSTPS 

PSTRS 

MSTPS 

PSTRS 

PSTRS 

PSTRS 

PSTRS 

PSTRS 

PSTRS 

PSTRS 

RSTPS 

RSTRS 

PSTRS 

PSTRS 

PSTRS 

RSTRS 

PSTRS 

HSTHS 

PSTRS 

RSTRS 

PSTRS 

PSTRS 

PSTPS 


ARl  FRC 
AHL  FRC 
APL  FRC 
ARl.  FRC 

ARL  FRC 

API  FPC 

ARL  FPC 

ARl  FPC 
API  FRC 
AFL  FRC 
ARL  FRC 
ARl  FRC 
AHL  FRC 
ARl  FRC 
API  FRC 
APL  FRC 
ARL  FRC 
ARl  FRC 
ARl  FRC 
ARl  FRC 

AF l  FRC 
ARl  rUC 
ARl  FRC 


nr>n  •oooonnoo 


Z2?  *  KP??K*(3.U«  ZIM-7.0*  Z(N*n*5.(l*  Z(N*2>-  7(N*3>> 
or?2  *  *.022**  n.o»cfi <m -7.o*dr  in* n  .s.o*dmn*?i -dr  t n ♦  a )  > 

C72?  *  M?2N*<3.(>»nZ(M-7.n*ltZ(N*l)*5.(i*D7(N*2)-l)Z(N*3)  ) 

GOTO  ? o 

If*  IF  (  IPCE’2  .NF.  4  .OR.  N  .IT.  N?R-|)buT0  15 
C  I«CF?a4  AM'  NaN2R-l  (M/2).  USF  NON-CFNTRAL  4  POINT  2ND  C I FFERFNCF  S 
•P  RESULT. 3*. PFSl  LT.4P 

IF(N  ,F(J.  N2R  .  ANC .  If*Cfc2  .fcO.  4 )  (i  fl  T  0  13 
IF(N  .Nfc.  Nlh  .OH.  IPCfl  .NF.  4 ) 0  0  T  0  1 A 
C  NaN 1 p  ANr  IflCE 1*4.  APPLIEO  FORCE  FN01. 

CALI.  FNCFPC 
AMS (NIP)  •  EMI 
GOTO  IP 

C  N  a  N2P  AND  IPCE?«4.  APPLIED  FOHCF  EN02. 

13  IF  *  If*CE  1  .ME  •  4 )  CALL  ENOFHC 
AMS(NPR)  *  £M2 
GOTO  IP 

•  P  MOTION. 11. MOTION. 18 

C  APPLIED  FORCE  (ANO  MOMENT)  AS  FND  CONDITION  4.  4/12/7* 

IFIIRCE1  .NE.  4  .OR.  N  .GT.  N1B*1)G0T0  13 
IF  (N  .GT .  Nit* )  GOTO  11 

c  NaNi  1 H  ANO  IHCE 1*4.  USE  01FFFRFNCE  IF (N* 1/2) -F (N ) ) / 1 OE TA?/?> . 

VR  a  2.n*PTC2M*(VS(N)«SNRM(N>.QS(K)*SNKM(N)-EFRl) 

V 7  a  2.P*RTL'2«*  ( VS  (N)  »SNKM(N)  *QS  (N)  *SNRM  (N) -EF71 ) 

GOTO  3S 
11  MaN-1 

C  N.MH.1  ANO  INCE2a4.  USF  0 IFFFRENCE  IF  (N*  1  /? > -E  IN- 1  /2 )  ) /PE TA? 

VRaRTn2M*(VS(M*SNPM(N)-0S(N)*SNKM(N)-VStM)*SNHM(M)*0S(M)*SNKMIM> 
V7aHTD?M* (VS  IN) »SNKM (N) *OS(N) *SNRM ( N ) - VS  I M ) *SNKM I M ) -OS ( M) »SNRM(M) 
GOTO  35 

13  IF  I IHCE?  .NE.  4  .OR.  N  .LT.  N2R-1IG0TC  31 
IFIN  .LT.  N2P ) GOTO  16 

C  NaN?p  ANr  IPCE2«4.  USE  DIFFERENCE  I F { N ) -E I N- 1 /2 ) ) / I Ofc T A?/? 1 

MaN?P— 1 

VR  a  2.0«RTC2M*(EFB2-VS(M)*SNRMIM)*OS(M)*SNKMiM) ) 

VZ  «  7.P*RTP7M*(EF72-VS(M)*SNKM(M)-QS(M)*SNF.M(.-'' ! 

GOTO  35 
16  M  *  M2H-2 

c  N»N?P-1  ANO  IHCF2a4.  USE  OIFFFPENCE  (F (Nal /2>-F (N-l/2) ) /PETA2 

VPaPTn?M*(VS(N)«SNRM(M-OS(N)*SNKM(N)-VS(M)*SNRM(M) *0S (M) *SNKM (M) 
V7aPTP2M* (VS (M •SNNM(N) «GS  <N> *SNHM ( N ) -VS ( M ) *SNKM ( M 1 -US (M) *SNRM ( M ) 
GOTO  35 


*T 

POUNOR.R 

*1 

IFMPCF1 
HOLNOII.  7 

.EH. 

4 )  GOTO 

20 

OAF 

IF ( I FCE 1 

.EC. 

4)  GOTO 

20 

•DFCN  ENDFRC 

SOBHPIIT  INF  FNOFHC 

THE  APL  FPC  CHANGE  TO  (1-0)  PFPSIL  FOR  R t AMS  HFQUIRFS  A  USER 
SUHPCUT INE  ENOFHC  TO  SUPPLY  FORCES  AND  MOMENTS  AT  NOTH  EM'S . 
FOR  AN  INITIALLY  STRAIGHT  REAM.  RIMaO.O.  7 ( N ) « ( N-2 ) OE T A? . 
FF71 >0  OEIRFASES  07(2).  FFHIXt  Dt CREASE S  DR(2) 

EF72>n  INCREASES  DZ(N28).  EFR2X)  INCREASES  DRIN2B) 

F. M 1  > 0  DECREASES  OR  ( 2 1  AND  INCREASES  DR(3) 

EM2>0  INCREASES  DR (N2H-1 )  AND  DECREASES  DRIN2B) 

CALL  MAIN 

TABULAR  INPUT  FOR  FF71.  EFTH.  ANO  EFMU.  10/27/7N 
MODIFIED  (12/7/76)  FORCE  RFL  ATIVE  TO  NORMAL 
FFR 1  a  EFMI:«EFZ1.  EM  «  7U*  (EF 7  1  *SIM EFTH )  *F  FR 1  <*COS  (EFTH)  ) 
r* ;  MR.'S  T  ON  L  A°tL  V  <  71  ,  ”F  7T  I  4("  .  TMF  7  («.(>)  ,fETi*t  uni  .  TMTN  <  40  I  . 

1  EFMUT(40) .TMMiKAO) 


AH  FRC 
AH  FRC 
ARL  FRC 
AH  FPC 
AH  FRC 
ARI  FhC 

APL  FPC 
APL  FRC 
API  FRC 
API  FPC 
API  FHC 
AH  FRC 
APL  FRC 
APL  FRC 
APL  FRC 
API  PRC 

APL  PRC 
APL  FRC 
APL  FRC 
AH  FRC 
AH  FRC 
API  FRC 
APL  FRC 
ARI  FRC 
APL  FRC 

arlfrc 

arlfrc 

APL  FPC 
API  FRC 
ARC  FRC 
•  L  FPC 
I  FRC 
mF  L  FPC 
AH  FHC 
API  FRC 
APL  FRC 
APL  FPC 
APLFRC 
ARLFhC 
API  FRC 

APL  FRC 

API  FPC 


FNDFRC 
FNDFRC 
f NCFRC 
ENDFRC 
FNCFPC 
ENDFRC 
FNDFRC 
F  N  OF  RC 
ENDFRC 

1P/2T/T6 

12/7/76 

K/27/76 

’/O/’T 

3/8/77 


U  U  »  »  <_>  O 


DATA  IFM>FR/o/ 

IFIIENOFR  .  GT  .  0  )  GOTO  1 
IENOFB  =  1 

BFAP  <*.  100) NF7PTS . (LAHtLV ( !  !  ♦ 1  =  1 .7) 

WRITE  (6.101  )  (LAHFLV  (  I  >  .  1  =  1 .7)  .NFZt-TS 

rf  An  efzt  <  ii  .tmfz < n . i  =  i .nfzpts> 

WR ITF  ( 6  •  1 03  >  (I. EFZT  ( 1 1  .  TMFZ  1 1)  .  I  =  1  .NF7PTS  > 
BFAO  (5.1001NFTHS  •  (LARELV ( 1 1 < I  =  ] <7 ) 

WPITF (6.104)  (LAHFLV  (  1  )  .1  =  1.7)  .NFTHS 
BEAD  (b.lOZH  EFTRTI  I) .TMTH< I) ,1  =  1 .NFTHS) 
WRITF (6.103) (I.EFTHT(I) .TMTHII) .1=1. NFTHS) 
BFAP  (5.1001NFMUS  .  (LARELV (11.1  =  1.7) 

WRITE <6. 105)  (LABEL V ( I ) .1*1.7) .NFMUS 
READ  (5.10?) (  EFMUT(I) .TMMU(I) .1*1. NFMUS) 
WRITE (6. 103)  (I.EFMIJT(I)  ,TMMU(T>  .1*1. NFMUS) 


1  CALL  nv(>INT(TIME.EF71.TMFZ.FFZT.NFZPTS.2> 

CALL  PVniNT(TIME.EFTH.T*TH,EFTHT.I'FTHS.2> 

CALL  PVPINT (TIHE.EFMU.TMRU.EFMUT.KFR05.21 
FFP1  =  FFPU*EFZ1 

F“l  =  -7U*(SIN(EFTH)»EF71  ♦  ( 1 . O-COS ( EFT*> ) *EFR 1 ) 

COPF  CHECKING  PtOTS.  10/2H/76 
P (53)  =  EF71 

BFP1 ACE  COMPUTATION  OF  FFZ1  ANP  EFZ?  IN  ENOFRC  12/7/76  CAPPS 

interpolate  for  normal  at  end  with  qijaopatic  interpolation 

XM*1  *  DE  -  PETA2 
XMX2  *  PE  -  2.0«0FTA2 

x-xo  =  PE.  X2-X1  =  Xl-xo  «  0ETA2. 

SNHFNP«  (XMX2*  (XMX1«SNR  (MB) -?.0*OF*SNB  «MB*1)  )  *DE*XMX1«SNR  (NIP*?) 

1  )*RP22M 

SNKFNPs (XMX2» (XMX1*SNK (N1B)-2.0«OF*SNK (M1P*1 ) ) »DE*X  MX  1 *SNK (NIP.?) 
1  )*PP?2M 

AAA  =  SORT (SNPFNO**2  ♦  SNKENO***?) 

SNRFNP  =  SNPENP/A A A 

SNKFNp  =  SNKENO/AA A 

TFFP1  s  EFR1*SNREN0  -  EFZl*SNKENO 

EFZ1  =FF R 1 •SNN  FNO  ♦  EFZ1*SNR£NP 

FFR1  3  TFFR1 

TFFP1  =  FFR1«SNH(MH)  -  EFZ1*SNK(N1H) 

FFZ1  =  FF01«SNK(Mfi)  ♦  EF  7 1 *SNR  (NIP) 

FFP1  =  TFFP 1 
PFTUPN 

)PC  F OPM A T(I10.7A10) 

101  FORMAT! • 1SUFR0UTINE  ENOFRC  10/27/76.  TABULAR  DATA  FFZl.EFTM.FFMU* 

1  /•  TAPI.E  1  .' .7A10.I5.*  P0!NTS*/*X.*I*.6X»*FF2T  *,9X.*TIME* 

102  F ORMaT ( ?E 1 5 • 7 ) 

103  FORMAT (IS. IP ?F 15. 7) 

]0*  FORMAT!/*  TARLF  2  .  *  .  7A 1 0 .  IS . •  RO I NTS • /* X •  •  I  * . NX . • F F THT *  .  9 X .  •  T  IME  • 
105  FORMAT  (  /  •  TAHIR.  3  .  *  .  7  A  1  0 . 15  .  •  POINTS*/**. »I».6X.*EF  MOT  **9X.*TIME* 
ENO 

’  OFNT  FAPFPC 


RPSLH’.a 

FXTFRNAL  WORK  FOR  APL  FHC  OF  4/12/70. 
MAIN.2P 

COMMON  OWFlP,|.ThlP.rwF?P.OTH?P 
RPSL 1 0. *3 
0TH1P  s  0.0 
0WF1P  *  P.0 
CWF2P  s  0.0 
OTH/o  =  0.0 

MOTION. *0 


(6/3/76) 


1  C/27/76 
10/27/76 
tf/27/76 
10/27/76 
10/29/76 
1C/27/76 
l C/27/76 
10/27/76 
1  C/29/76 
10/27/76 
10/27/76 
10/27/76 
10/29/76 
10/27/76 
10/27/76 
10/27/76 
10/27/76 
10/27/76 
10/27/76 
10/27/76 
1C/27/76 
10/27/76 
10/29/76 
CC/10/76 
7/20 
FRCP21 
E  R  0  U  2 1 
FRC021 
FRC021 
FROD21 
FROP21 
FROD21 
FR0D21 
FR0C21 
FRU021 
FR0021 
ER0U21 
ER0021 
F  R U021 
EROC21 
12/7/76 
12/7/76 
12/7/76 
1 0/27/76 
1  C/27/76 
.10/27/76 
) 10/27/76 
10/27/76 
10/27/76 
)  1 1'/?7/76 
) 10/27/76 
10/27/76 


FAP  FRC 

FAP  FRC 

FAP  FPC 
FAR  FRC 
EAR  FRC 
FAR  FRC 


124 


•  n  r>  »  •  o  •  o  *  o 


IF (N  .GT.  Nlfl)GOTO  130 
0PI1  »  CR1 
DZI1  *  071 

130  CONTINUt 

IFIIPCF1  .NF.  4)  ROTO  131 
NaNlR 

OwF  *-0.5*(EFHl*(r«Il*OP(M )«EF71»tl>ZIl*0ZIN> ) ) 

R2*RT02*  <-3.0*<R IN) «0R IN) ) *4.0*  |H (N.l ) »OH <N*1 ) ) - (P IN*2)  *DR IN*?) >  > 
Z?aPTU2* 1-3.0* IZIN ) ♦OZIM  >  *4 . 0* ( Z  IN *  1 ) *02 IN*1 > ) -(2 IN*2) *DZ IN*?) > > 
0  *  SORT ( 72**2  ♦  H?**?> 

SNRN  a  72/0 
SNKN  «  -R2/0 

OTHFT  a  ARCS IN I  SNK IN  1 *SNRN  -  SNR (N ) *SNKN) 

OWN  a-o.S*FWl*(OTMlP  ♦  OTHET) 

TNPO  a  TNRS  ♦  0.5MOWF  ♦  OWM  ♦  DWE1R) 

DNE1P  a  OWF  ♦  OWN 
DTH1P  a  OTHET 

131  IF  I  IPCF.2  .NF.  4 )  GOTO  132 
NaN?R 

DWF  a  0.5«<EFH2*I0RI.0RS)  ♦  EF 72* IDZ I *OZS )> 

R2aPTD2* IR IN-2) .OR  IN-2 1-4.0* (R(N-1 ) ♦OH(N-l) ) *3 . 0* CP  I N > ♦OH (N ) ) ) 
Z?aPT0?a(Z(N-2) *02  IN-2 ) -4 . 0* I Z I N- 1 ) ♦OZ(N-l) I ♦ 3. 0* ( Z (N ) *02 (N) > ) 

I)  a  SOP T 1 22**2  *  R2**2I 
SNRN  a  72/C 
SNKN  a  -R2/D 

OTHFT  a  ARCS  I N I SNK < N ) *SNRN  -  SNP|N)*SNKN) 

OWM  *-0.5*EM2*lOTH2P  ♦  OTHFT) 

TNRG  a  TNRG  ♦  0.5*<OWF  ♦  Own  ♦  OwF2P> 

0WE2P  a  OWF  ♦  OWN 
0TH7P  a  OTHET 
13?  CONTINUE 
•IOFNT  11/2/76 
*1  RFSL1D.5 

INSERT  TO  REMOVE  FAILURE  HALT  AT  FREE  OR  FORCED  FNO •  11/2/76 

I  RMSTRS.13 

INSERT  TO  REMOVE  FAILURE  HALT  AT  FPFE  OR  FORCED  END.  11/2/76 
IF  <02?  .GT.  0. 0 ) GOTO  2 
IF ( IP  .EQ.  2) GOTO  2 

IFIN  .EQ.  N1R  .AND.  I HCE 1 .EQ.4 ) G??a 1 . OE-1 0 
IF  IN  .FQ.  N?R  .ANO.  IBCE2 .F0.4) G??» 1 . OE- 1 0 
2  CONTINUE 
T  OF  NT  INVFLC 
I  INVEL.6 

LwTFPAL  VELOCITY  .  VH*  TONARO  FNU  1. 

ENTIRE  ROO  MOVING  KITH  CONSTANT  VELOCITY.  7/20/76 
WRITE  16.333) 

333  FOPHATI//'  LATERAL  VEIOCITY  TOWARD  END  1.  l/P/76*//) 

UO  29  NaN 1 V  «N2V 
DZ(N)a-VR 

29  WRITE  16,666) N.OZ <M .DR (N) 

RETURN 

666  FORMAT!*  Na*,I3,*  OZ(N)a  *.1PF13.6.'  ORIN)a  *,F13.6> 

IDENT  ERO026 
I  RPSL10.S 

EROSION  OF  FNO 1  •  12/26/76.  PUT  AFTER  APLFRC  AND  SHR3/1. 

INCLUOE  SUHHOUT I NF  fcROOE  ANO  CORRECT  SUPROUTlNF  ENOFWC. 

I  MAIN.2P 

COMMON  OE 
•I  PPSL 1 0. 1 ?Q 

CALL  EROOE 
•I  RESULT. 61 

IFIN  .Fo.  NIP) VS(N) a ( AMS IN.l ) -AMS  IN) ) / ISRA* <0ETA2-fF) * 

IFIN  ,  KC.  N 1 H  I  (jOTC  200 


EAR  FRC 
EAP  FRC 
EAP  FRC 
FAR  FRC 
EAP  FRC 
EAF  F PC 
EAP  PRC 
EAP  FPC 
EAP  FRC 
EAP  FRC 
EAP  FRC 
EAP  FRC 
EAP  FRC 
EAP  FRC 
EAP  FRC 
EAP  FRC 
EAP  FRC 
EAP  FRC 
EAP  FRC 
EAP  FRC 
EAP  FRC 
EAP  FPC 
EAP  FRC 
EAP  FPC 
EAP  FRC 
EAP  FRC 
FAP  FRC 
EAP  FRC 
EAP  FRC 
FAP  FRC 
EAP  FPC 


11/2/76 

11/2/76 

11/2/76 

11/2/76 

11/2/76 

11/2/76 

11/2/76 


l/H/76 

1/8/76 

1/8/76 

7/20/76 

1 /P/76 

1/8/76 

1/8/76 

1/8/76 


F  POD 


F  F  CD?  1 
EP0021 


C'  u  u  u  u 


*1  MOTION. 26 

IF(fc  ,N£.  Nit*  .OR.  IRCF1  .NF.  4100TU36 

C  REPLACES  VP.VZ  f  OMPUTAT  I  ON  IN  NOTION  aHFN  NaNlP  ANC  IHCE1*4 

VU  s  (  VS  <IM)  «St  RM  (N  )  -OS  (N)  «SN*  «  (N)  -F  FP  1)  /  (  0  ,S«fifc  TA?  -DM 

V7  a  ( vs  (N)  *SNKM  (N  )  »OS  <N)  *SNRM  (N)  -FFZ1 )  /  (  O.S«DETA?  -  lit) 

36  CONTINUF 
•I  STRAIN. fll 

I F  ( Jl  .GE  .  MR)  GOTO  44 
FPSS1 (I)  a  FPSS2 ( I )  a  n.O 
FPSANfi(t)  a  EPSANR(I)  a  -1.0E20 
GOTO  46 
44  CONTINUE 
•I  STRAIN. V7 

IFINOl  .bE .  MR)  GOTO  4B 
01  a  02  a  o.O 
GOTO ( 71 .50) .L INK 
4R  CONTINUE 
*1  POATA.7 

MR  a  2 
*0  POATA.M 
100  Nlfl  a  N 1 A 
RFTUPN 

*  AF 

•OF  OK  ERODE 

SURPOUTINF  FPOCfc 

INSFPT  ANY  INPUT  FOR  ERODE  AFTFP  INPUT  FOR  ENDFRC • 

THIS  USER  SUBROUTINE  COMPUTES  THE  FROSION  ON  FNU1  OF  A  ROD. 
THIS  VERSION  ASSUMES  TABLE  OF  FROSION  FROM  FND  OF  ROD  VS  TIMF. 
TRF  TABLE  PPESUHABLY  AGREES  N I TR  ENDFRC  VS  TINE  TABLE  IN  ENDFRC 
AND  STRESS-STRAIN  TABLE  READ  IN  INPUT. 

•CALL  MAIN 

DIMFNSION  LAHLTE (7) .TEROO (100) .  XEFO ( 100) 

OATa  ItPOOE/O/ 

IF ( IEPOOE  .Efl.  1 ) GOTO  10 
OECON  a  DETA2/4 .0 
I  ERODE  a  \ 

FROOP  a  0.0 
OF  a  0.0 

READ (5. 1001NTFPTS. (LAPLTE(I) .la], 7) 

NPITE (6.101 INTFPTS. (L ARl TE < 1 1 . I  a \ ,7 > 

PFAO  (5.102) (  TFROO(I) .XENDII) ,I*1.NTFPTS) 

WRITE (6, 103)  (  l.TFROO(I) ,XFNO< II . I«1 .NTFPT5) 

100  FORMAT ( I  1 0 , 7A 1 n ) 

101  FOPHAT ( *  1 SUPPOUT INF  FPOOE  5/26/77*  TABULAR  OATA  TINE  VS  XEND*/ 

2  IS.*  POINTS*. 5X.7Alo,/4X.*I«,7X.»TINF».Rx.»XEND») 

102  FORMAT (?t 15.7 ) 

103  FORMAT) I5.2F15.7) 

10  CALL  OVD  I  NT  (TIME  •  EROO  .TER()O.XEND.NTF.PTS.2> 

DFLOF  a  FPOn  -  ERCTP 
IE  a  OF  ♦  OFLOF 
F.POOP  a  EROO 
P(54)  a  F POD 

C  COOF  CHECKING  PLOT  5/27/77 

1°  IF (OF  .LT.  OFCOM (GOTO  ? 0 
OF  a  OF  -OFTA2 
P(N1B)  a  Z (N 1 R )  a  -1.OE20 

C  THE  FOLLOWING  CARO  IS.NEFOEO  IF  ThF  DECK  INCLUOFS  MSOSVS. 

AMS(NIH)  a  IISIMHI  a  VS(NIB)  a  fl.p 
N1R  a  N 1  v  a  N1A  a  MR»1 
GOTO  )R 
20  PF  TURN 
F  NO 


FP0C21 

tFor.2i 

FPOD21 


FF0D21 

M-0021 

EF0D21 

FP0D21 

FR0D21 

FF0D21 

FP0D21 

FP0D21 

FP0021 

FP0D21 

FP0D21 

6POD21 


FRCD21 
FF0C21 
EF0C21 
FR0026 
FP0D26 
t  P0026 

t  R0D26 
FRO021 
FPPD2I 

FR0D21 

FP0D26 

EP0D26 
FP0D26 
f  R0D26 
FFUD26 
FPO026 
EP0D26 
FRCD26 
FR0026 
FFCD26 
FBUD26 
FP0U26 
t  R0D26 
FFCC26 
C  C 
CC 

FPCD21 
FPCD21 
FPUD21 
FR0D21 
FP0D21 
FPO021 
FI-0021 
FP0U21 
F  F  0021 


126 


HFPfcAT  OF  6/P/77  PUN 


FPOD1M3  BOD  <254  PY  8.126  MM) 
36  b  2.0 

130  o  ooou  .snonnoF-b 

4  4 

n  0  booq  60.44757  .001 


1 .50571 4E  1 1 

.27 

1.390( 

0.007 

1  .390 

F  9 

0.00A 

1.510 

09 

0.009 

1.564 

F9 

0.010 

1  .597 

E9 

0.011 

1.624 

f  9 

0.012 

1  .640 

F9 

0.013 

1  .669 

F9 

0.014 

1  .680 

F9 

0.015 

1.705 

£9 

0.016 

1.721 

1 9 

0.017 

1  .736 

£9 

0.010 

1.750 

E9 

0.019 

1.763 

F  9 

0.020 

1.776 

t  9 

0.022 

1  .799 

F9 

0.024 

1 .020 

F  9 

0.026 

1.040 

£9 

0.020 

1  .060 

F9 

0.030 

1  .674 
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O  6  10  16  20  P5  K  SS  MO  MS  SO  66  60  66  70 

TINE  IHICF06ECONDS> 

urmiwQ  ETBB  .  .Q40  ,  7.6331  lajrE* 

_ _  ETfU  COMPONENT 

_ ETPE  COMPONENT 

Figure  E.2.2.  Strain  plots  at  19.26  mm  and  39.67  mm. 
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APPENDIX  F 


SOME  ALTERNATE  USER  SUBROUTINES 


UPDATE  input  cards  for  alternate  versions  of  the  user  subroutines  INGEOM, 
PRESS,  and  ENDFRC  are  listed  in  this  appendix.  Other  versions  of  the  user 
subroutines  are  listed  in  Appendices  C  and  E.  Descriptions  for  user  sub¬ 
routines  are  given  in  Section  4. 

The  following  are  listed  in  the  order  given: 

PRESS  for  Two  Phase  Pressure  Decay  in  a  Cylinder.  This  is 
described  in  Section  4.2.3. 

INGEOM  for  an  axisymmetric  cylinder.  This  is  described  in 
Section  4.1.2. 

ENDFRC  for  one  or  both  ends  free .  This  is  described  in 
Section  4.5.1 . 
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r  I  OFF  SUN  7SI  (50)  .  TSI  (5  0)  .RSI  1 50 > • TS 1 1 03 ) ,PS 1 1 03) . TSP 1 1 03 ) • 

1  PSP( 103)  .PTC  (p) 

TATA  IPpESS/O/ 

IFflPHfSS  .6T.  It )  GOTO  30 
IPPFSS  c I 

GFAP (5.  10  0)  PPAH.pETA. (P10< I) .1*1. t) 

PFA0(5.105)CCNCZ»CCNCT«C0NCP.TDIF 

VPITFIlS.  10?)  PPAR.flETA.  <PI0<  I)  ,  I*l.f>)  .TDIF.C0NC2.C0NCT  .CONCP 
»FA0(5. 103) ALPHAS. POOSE .ISP  12/2/75 

PE*PI5.101) IZST II) *TSI 1 1) .PSI II) .lal.ISM) 

RPITEtE>.tO*)  ALPHAS, I7SIII) .TSI  II) .PSI (I) .1*1  .ISM)  12/2/75 

C  CONVERT  CATA  TO  Z*IN.  T»SFC.  PaLts/SCIN.  12/2/75 

PC  2  1*1  ,  ISP 
ZSI  (I)  «  C0KC7*ZS 1(1) 

TSI(I)  *  CONCT*(TSI (I)-TOIF) 

PSI (I)*CCNCP«PSI (1)  12/2/75 

?  CONTINUE 

WPITF (A. 104) (ZSI (!) .TST (I) .PSI II) »!»1.ISM) 

C  IMTIATE  DATA  FCD  CCMPUTlMi  THF  PPESSUPE 

CO  ?«  N*MV.N2V 
DEN  m  FLOAT (K-?)*CETA? 

CALI  POO! NT (CFN.TS(N) .2SI.TSI.ISP.2) 
call  ovr  jntiofn.ps ipi)  .zsi.psi.ism.?> 

TSHIN)  *  TSIN)  -  ALOUIPHAR/PSIN))  /  ALPHAS 
PSIM  *  PSIN)*EXP(AlPpaS«TSIN)) 

PSP(N)  *  Po 4h »r 4p (PfT A*TSP ( M ) 

?H  CONTINUE 

VRITF (H. 1 1 1 ) IZIN) ,TS IP) «PSCN) .TSRIK ) .PSB  IN) ,N*N 1V.N2V) 

COMPUTE  Thfc  PRESSURE 

30  £  ALPHT  a  FXP (-ALPPAS«T IMt I 
FhETAT  *  F*P(-R£TA«TIPF) 

CO  3«  HM  V.K?V 
IPITIMF  .LT.  TSINHOOTO  40 
IFITIPF  ,(5T  •  TSH  (N)  )  ROTO  31 
PIN)  »  -FALPHI*PS(M 
FOTO  30 

31  PIN)  *  -EHFTaT*pSP IN) 

30  CONTINUE 
40  RETURN 

1 00  E  OOMA  T (cK 10.3.NA1P) 

101  EUHPATI3M0.3) 

1 Oc  F0PV4TI4E 10.3) 

102  FOPNATI//  IX'PRESS  SUPPOUTINE  FOP  1-P  STRUCTL HE . 3/ 1 0/7F . • 

1  /.  FPCM  TAPLFS  OF  ARRIVAL  TI*«E  AND  PEAK  PRESSURE  VS  POSITION.  •/ 

?  •  POSITION  a  0  f CROFSP0N0S  TO  N*2  UP  1-0  REPSIL.V 

3  •  THF  (.IASI-STATIC  PRESSURE  I S » 1 PF  1 5 . 7 .5 A .  •  w I Th  EXPONENTIAL  » 

4  •  tAMPING  COEFFICIENT  *«F15.7/10X.6A10/ 

*  *  INPUT  TIPF.  TSl .  OtPUCfcD  PT  TIME*'  .F15.7/ 

5  «  COPVFRSIOK  COEFFICIENTS  FOR  Z.  T.  AND  P  ARE  »3E15.7) 

103  FORMAT ( ?E 10.3,110) 

104  FORMAT I//20X, *  INPUT  TATA  •/•  ALPHAS* « , 1  PE  1 5 . 7 . 

t  //7X,  *7SI  *.!?*•  'TSP.12X.  »PSI'/(3E1S.7)  ) 

10h  FORMAT 1//SX, *TAMULAO  taTA  AFTER  ("INVERSION*/ 

1  7X.*7ST* .12X.* TSI • ,1?X . »PSI*/( 1PE3E 15.7) ) 

111  FORMATI//SX. FINAL  IATA  » /H « * Z • 1 4 X • TS *  1 3X *PS • 1 3X 

1  'TSP* 1?> ‘PSP*/// ( IPSE  15.7) ) 


*P  IK0FUP.7.  INGECP.aT 
*r 4i  l  ►  «fp 
r 

r  ikpi-.jm  -ci-  air'vn'tT^r  rvL  'kiftrk .  ';•/  7/75. 


c 

C  T  PHOGhAp  Til  USE  GAUSSIAN  IK TF  uH AT  t ON 

TGAt'SS  a  1 

RFAIHN.IOOICYLL, RADIUS 

w*ITF(f ,101)CYLL*PAim.S 
0FTA1  a  2.0A?.1MS<52N535PR7P3*PADII  S 
C  0FTA1  «  2*PI*RA0IUS  GIVES  ENFPGY  FOR  FNTIRE  CYLINOFP 

l)FT A?  a  CYLL/FLOAT(NPRSH) 

If)  ■  0 

C  *****  F VALLATE  PIN)  AKO  2  (N  )  ♦  K«MP.N2fc«  ***** 


DO  10  P»?.N?S 

R(N)a«ACIUS 

2  (N) aFLOAT (N-2) *OFTA? 

10  CONTINUE 
RETURN 

100  FOPPAT (2E12.F) 

101  FORpAT  ( //  •  TNGF.OP  FOR  CYLINOFR  12/3/75.  LENGTH  ■  •  1 PF 1 5 .6  *5X  . 

1  •  RAO IOS  a«,FlS.A/) 

EM) 


*AF 

*OFCK  FMPFPr 

SI'OROI'TInF  FNDFPC 
•CALL  PAIN 

THF  APL  FPC  CHAM5F  TO  (1-0)  PEPSlL  FOR  HEAPS  RFCUIRES  A  USFP 
SUPRC’UTTNP  FNOFBC  TO  SUPPLY  FORCES  AND  POPENTS  AT  HOTH  FFDS. 
FOP  AN  INITIALLY  STRAIGHT  REAP,  R(N)a0.0«  2  ( N )  ■  ( N-2 )  OE  T  A?  , 
FF21 >n  ofcheasFS  02(2),  EF«I>n  OFCRFASFS  DRI2) 

FF7?>n  IKCPFASFS  07IN2H),  EFP?>C  INCREASES  DP(N?P) 

F  P I  >  (I  OFCPFASES  np(2)  ANO  INCREASES  DR(3> 

F‘«?>n  INCREASES  0R<*2*-1)  ANP  PFCHEASFS  OP(N2H) 

THIS  VFRSTON  IS  FOP  A  RFAP  PITH  ON F «  TWO,  OR  NO  FRFF  FNOS, 
TATA  TFK'FO/n/ 

IF ( tEKPFP  ,GT.  OIPETUPN 
IFMFR  «  1 
FFO)  a  n.il 
FF71  a  0.0 
API  a  O.n 

FFU?  a  O.n 
FF2?  a  n.n 
F«?  a  n.o 
PFTURF 
FK'O 


FFOFRC 

ENTFRC 
FNDFRC 
ENDFRC 
FNDFPC 
ENDFRC 
LKCFhC 
FI TFPC 
ENT FPL 
f KCFHC 
FM'FPC 
FI  DP  PC 
EIOFRC 
EIOFPC 
EICFRC 
FNPFRC 
EIOFRC 
ENDERC 
FNDFRC 
EIOFRC 
FNDFRC 
EIOFRC 


( 
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STBFM 

»p«l  in.s 

'F'-  ttT°4rN  FnFCC.V  CONH'lTAT  If,H  Fi’H  A  iJEAM,  1/20/74  STHEN 

t^'ST-F.T 

Nf-  STRAIN  t-NHHPY  COMPUTATION  1/20/76  STHEN 

PImFNSIUN  SST  <5>  .0FST(«>)  .LST  (S)  STHEN 

h*>STRS.33 

npm  STRAIN  FNFPGY  COMPUTATION  1/20/76  STktN 

I F  ( IP  .Ft..  2  )  OOTii  t»03  STHEN 

SST(J)  r  STR2?  STHEN 

LST ( J)  a  J  STHEN 

PMSTPS.36 

NEW  STRAIN  FNFCC-V  COMPUTATION  1/20/76  STHEN 

HMSTRS.3H 

NFW  STRAIN  FNFWPY  COMPUTATION  1/20/76  STPEN 

IF  ( IP  ,FC.  2 )  fiOT 0  ?S  STHEN 

IFINSFL  .GT.  1 ) GOTO  10  STHEN 

OSTPFM  * <  SS22 I **2  STHEN 

GOTO  3"  STHEN 

10  IF ( SS22  .UE.  0 . 0 ) GOTO  12  STHEN 

S<522  a  -SS22  STPEN 

00  11  l.»)  .NSFL  STHEN 

11  SST  ( 1. 1  a  -SST(L)  STHEN 

1?  00  13  l. *1  .NSFL  STHEN 

13  CFSTU)  a  (SI GZZtl.)  -  SST(L))/F  STHEN 

LLASTaNSFL-1  STHEN 

TO  IS  L=1,LLAST  STHEN 

JSTAPT  a  U»1  STHEN 

00  14  JajSTART.NSFL  STREN 

IFtCFCTtLl  .IT.  nEST( J! 1GPT0  U  STHEN 

OSTFMP  a  OEST(L)  STHEN 

OFSTtLl  *  OEST(J)  STHEN 

rFST(J)  a  OSTF**P  STHEN 

ISTFVP  a  LST(L)  STHEN 

LST(U  a  LSTIJl  STHEN 

IST(J)  a  LSTFMP  STHEN 

14  CONTINUE  STHEN 

15  CONTINUE  STHEN 

Fm  *  1 .n  STHEN 

TSThfn  a  n  .  (i  STHEN 

00  ?c  Lal.NSFL  STHEN 

FSTAH  a  F««*E  STHEN 

OFPSTP  a  -SS22/FSTAH  STHtN 

OSTRFn  a  nsT-fcN  ♦  <SS22)»«2/F*«  STHEN 

IEIOFPSTH  .LE.  OEST(LI)  GOTO  30  STHEN 

COPNFP  AT  rFST(L)  If  UNLOAOINtt  Cl'HVE  STHEN 

SS22  a  -FSTAH*  (OE'PSTO-rEST  u. )  1  STHEN 

p.sTHFN  a  0STHE4  -  (S«?2I**2/F«  STHEN 

LI  a  LST (L)  STHEN 

E M  a  LN  -  HTILll  STHEN 

0FST1  a  CFST(l)  STHEN 

r<>  ?0  jal  .»>SFL  STPFN 

OFST(J)  a  IIFSTIJI  -  OFSTL  STHEN 

?0  CONTINUE  STHEN 

2?  CONTI NUF  STHFN 

30  STHFN  a  STUFF  ♦  S0G*M  f N I *05 THEN  STPEN 

3S  CONTINUE  STHEN 


LIST  OF  SYMBOLS 


[FORTRAN  name  in  brackets] 


r  position  vector  to  a  point  on  the  reference  surface 

R,Z  position  components,  r  =  R  i.  ^  +  Z  i_2  [R(N),Z(N)] 

t  time  [TIME] 

area  or  weight  associated  withc^  [W(K)  ] 

At  time  increment  [DELTAT] 

2  2 

A£,A£;  spacing  of  C  between  mesh  points  [DETA2] 

5  5  is  52 

12  12 
5  ,£  material  (Lagrangian)  coordinates  (called  n  ,n  in  Reference  2) 

S  normal  distance  frim  reference  surface 

C  at  k'th  integration  station  [ZETA(K)] 

)  t,  on  lower  (upper)  surface  of  shell  [ZL,(ZU)] 

■JC  U 

v  Poisson's  ratio  [FNU] 

p  mass  density  [RHO] 
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